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Magnetic  nanostructure  has  been  a  new  trend  because  of  its  application  in  mak- 
ing magnetic  sensors,  magnetic  memories,  and  magnetic  reading  heads  in  hard  disks 
drives.  Although  a  variety  of  nanostructures  have  been  realized  in  experiments  in 
recent  years  by  innovative  sample  growth  techniques,  the  theoretical  study  of  these 
devices  remain  a  challenge.  On  one  hand,  atomic  scale  modeling  is  often  required 
for  studying  the  magnetic  nanostructures;  on  the  other,  these  structiures  often 
have  a  dimension  on  the  order  of  one  micrometer,  which  makes  the  calculation 
numerically  intensive. 

In  this  work,  we  have  studied  the  electron  transport  theory  in  magnetic  nanos- 
tructiures,  with  special  attention  to  the  giant  magnetoresistance  (GMR)  structmre. 
We  have  developed  a  model  that  includes  the  details  of  the  the  band  structmre 
and  disorder,  both  of  which  are  both  important  in  obtaining  the  conductivity.  We 
have  also  developed  an  efficient  algorithm  to  compute  the  conductivity  in  magnetic 
nanostructures.  The  model  and  the  algorithm  are  general  and  can  be  applied  to 
complicated  structinres.  We  have  applied  the  theory  to  current-perpendicular-to- 
plane  GMR  structures  and  the  results  agree  with  experiments.  Finally,  we  have 
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searched  for  the  atomic  configuration  with  the  highest  GMR  using  the  simulated 
annealing  algorithm.  This  method  is  computationally  intensive  becaiise  we  have 
to  compute  the  GMR  for  10^  to  10^  configurations.  However  it  is  still  very  efficient 
because  the  number  of  steps  it  takes  to  find  the  maximum  is  much  smaller  than  the 
number  of  all  possible  GMR  structures.  We  found  that  ultra-thin  NiCu  superlat- 
tices  have  siurprisingly  large  GMR  even  at  the  moderate  disorder  in  experiments. 
This  finding  may  be  useful  in  improving  the  GMR  technology. 
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CHAPTER  1 
INTRODUCTION 


Magnetic  nanostructures  [1-2]  are  magnetic  devices  which  have  at  least  one 
dimension  of  the  order  a  nanometer,  10~^  m.  These  devices  have  apphcations  in 
making  magnetic  sensors,  magnetic  memories,  and  magnetic  reading  heads  in  hard 
disks  drives.  Predicting  the  electrical  transport  properties  theoretically  is  a  chal- 
lenge because  the  detailed  arrangement  of  atoms  affects  the  observable  properties 
significantly.  In  this  dissertation,  we  will  develop  methods  to  predict  the  transport 
properties  of  magnetic  nanostructm-es  by  taking  into  account  band  structiue  and 
disorder.  We  will  first  develop  models  to  compute  the  giant  magnetoresistance 
ratio  (OMR),  which  is  the  fractional  change  of  the  resistance  due  to  an  applied 
magnetic  field,  of  multilayer  structures  and  compare  our  calculations  with  exper- 
iments. We  will  then  solve  the  "inverse  problem,"  which  is  to  find  the  atomic 
configviration  that  has  the  highest  GMR. 

Although  the  method  we  develop  can  be  applied  to  various  geometries,  in 
this  dissertation  we  concentrate  on  magnetic  superlattices,  which  are  multilayers 
formed  by  repeatedly  stacking  magnetic  layers  followed  by  non-magnetic  layers. 
The  thickness  of  each  layer  is  of  the  order  of  one  nanometer;  however,  the  other 
two  dimensions  can  be  as  large  as  one  micrometer,  10"^  m.  This  geometry  can  be 
grown  by  techniques  used  in  depositing  thin  films.  In  the  following,  we  will  first 
discuss  the  GMR  experiments,  explain  several  pictures  for  understanding  the  gi- 


1 


2 


ant  magnetoresistance,  and  discuss  the  recent  trend  of  optimization  on  the  atomic 
level. 

1.1    GMR  Experiments  ' 

The  first  "giant"  magnetoresistance  measurement  was  that  of  Baibich  et  al.  [3], 
where  they  observe  an  almost  50%  reduction  in  the  resistance  of  Fe/Cr  superlattices 
in  an  applied  magnetic  field  of  20  kG.  As  shown  in  Figure  1.1,  the  resistance  at 
zero  field,  Rq,  is  the  highest.  As  the  magnetic  field  is  applied,  the  magnetization 
of  the  adjacent  Fe  layers  align,  which  causes  the  resistance  to  decrease.  At  a  field 
larger  than  the  saturation  field  Hs,  the  magnetization  of  the  adjacent  Fe  layers  is 
parallel  to  each  other.  Therefore,  the  resistance  stays  at  the  value  of  the  saturation 
field,  Rsax.  The  giant  magnetoresistance  ratio  is  defined  as 

GMR=  (1.1) 

In  this  experiment,  the  current  travels  in  the  plane  of  the  multilayers,  which  is 
referred  to  as  the  ciurrent-in-plane  or  CIP  giant  magnetoresistance. 

In  this  dissertation,  we  focus  on  another  geometry,  the  current-perpendicular- 
to-plane  (CPP)  geometry,  in  which  the  current  fiows  perpendicular  to  the  plane  of 
the  multilayer.  The  CPP-GMR  is  usually  larger  than  the  CIP-GMR  for  the  same 
lattice.  Most  experiments  to  date  are  for  the  CIP  geometry,  while  a  smaller  portion 
of  the  experiments  are  for  the  CPP  geometry  [4-15].  Current-perpendicular-to- 
plane  experiments  are  difficult  to  perform  because  the  thickness  of  a  superlattice 
is  of  the  order  of  a  hundred  nanometers  at  most,  while  the  area  of  the  layer  is 
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Magnetic  field  (kG) 


Figure  1.1:  Magnetoresistance  of  three  Fe/Cr  superlattice  at  4.2  K  [3].  The  current 
and  the  apphed  field  are  along  the  same  [110]  axis  in  the  plane  of  the  layers.  Hg 
is  the  saturation  field  required  to  align  the  magnetizations  of  Fe  layers.  As  the 
magnetizations  are  aligned,  the  resistance  is  reduced  from  the  peak  at  H  =  0. 

usually  much  larger  than  1  //m  x  1  //m.  When  a  current  passes  though  the  system 
perpendicular  to  the  plane,  the  resistance  is  too  small  compared  with  the  resistance 
of  the  lead  and  the  contact,  and  is  thus  hard  to  measure.  To  measiue  the  CPP 
resistivity,  one  must  improve  the  aspect  ratio  of  the  sample.  Small  cross-sectional 
area  to  length  ratio  has  been  produced  by  using  micro-structiured  pillars  (Figiire 
1.2)  [4],  by  growing  superlattices  inside  the  pores  in  polycarbonate  (Figure  1.3) 
[15],  and  by  connecting  several  superlattices  with  superconducting  leads  (Figure 
1.4)  [6-7]. 
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substrate 


ib) 


'  Mo    100  nm 
■  AI2O3  200  nm 
Mo    100  nm 
Au     300  nm 

Fe/Cr  500  nm 
Cr     *00  nm 

■Si02 


(c) 


Polyimide 


Figure  1.2:  Schematic  diagram  the  samples  growth  process  used  by  Gijs  et  al.  [4]. 
The  cross  sectional  area  S  ranges  from  6  to  130  /im^.  - 

1.2    Models  of  the  CPP-GMR  • 


There  has  been  important  progress  in  understanding  the  CPP-GMR  effect. 
Classical  resistors-in-series  models  can  be  applied  to  study  the  CPP-GMR  when 
the  layer  thickness  is  much  larger  than  the  Fermi  wavelength  [16].  The  origin  of  the 
spin  anisotropy,  which  measures  the  difference  between  the  resistivity  in  the  two 
spins,  is  more  difficult  to  understand.  Various  free-electron  models  [17-21],  tight- 
binding  models  [22-26],  and  density- function  methods  [27]  have  been  performed. 

In  the  following,  we  discuss  several  models  for  the  CPP-GMR.  In  the  CPP- 
GMR  system,  each  electron  has  to  go  through  every  layer  and  interface.  It  is  thus 
easier  to  understand  intuitively  than  CIP-GMR  system,  where  electrons  travels 
along  the  layers,  while  at  the  same  time  going  across  interfaces  of  layers  and/or 
scattering  from  the  interfaces.  ^ 
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Figure  1.3:  Schematic  representation  of  the  array  of  nanowires  in  the  insulating 
polymer  matrix  [15]. 

1.2.1    Resistors-in-Series  Model 

When  the  layer  thickness  is  smaller  than  the  spin-diffusion  length  [28-29],  one 
can  treat  the  two  spins  independently,  and  this  approach  is  called  the  two-fluid 
model  [16].  When  the  layer  thickness  is  much  larger  than  the  Fermi  wavelength, 
the  GMR  effect  can  be  explained  with  the  classical  resistors-in-series  model,  which 
assumes  electrons  pass  through  spin-dependent  resistors,  as  shown  in  Figure  1.5. 
Let  the  resistance  for  the  majority  and  minority  spin  in  the  bilayer  formed  by  a 
magnetic  layer  and  a  non-magnetic  layer  be  R  and  r,  respectively.  At  saturation 
field,  the  magnetic  moment  of  adjacent  magnetic  layers  are  parallel  to  each  other. 
Assuming  the  spin  diffusion  length  is  much  larger  than  the  layer  thickness,  the 
resistance  for  the  majority  and  minority  spin  is  R  +  R  and  r  +  r,  as  shown  in 


Nb 


I     I  SiQ, 


(a) 


(b) 

Figure  1.4:  CPP  sample  of  Cyrille  et  al.  [6].  (a)  Schematic  cross  section  of  a  CPP 
sample;  (b)  optical  micrograph  of  a  typical  sample. 

Figure  1.5(a).  Therefore,  the  net  resistance  for  the  saturation  field  is  given  by 


_    1  Rr 


2R  +  r 


(1.2) 


At  zero  field,  the  magnetic  moment  of  adjacent  magnetic  layers  are  assumed  to  be 
anti-parallel  to  each  other,  as  shown  in  Figm-e  1.5(b).  A  majority  spin  electron  of 
the  first  magnetic  layer  will  be  in  the  minority  spin  channel  of  the  second  magnetic 
layer.  Therefore,  the  resistance  for  each  spin  channel  is  r  +  i?.  As  a  result,  the  net 
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resistance  at  zero  field  is  given  by 


i?o  — 


R  +  r 


(1.3) 


Since  (i?  +  r)^  >  4Rr,  one  can  see  that  Rq  is  larger  than  i2sat,  which  is  the  explains 
the  GMR. 


(a) 


(b) 


R 


R 


R 


R 


Figure  1.5:  Simple  resistor  model  for  CPP-GMR.  The  circuits  corresponding  to  (a) 
parallel,  and  (b)  anti-parallel,  magnetic  configuration.  In  both  cases,  the  majority 
and  minority  spin  electrons  pass  through  the  left  and  right  series  of  resistors, 
respectively.  In  (a),  the  resistance  for  the  majority  and  minority  spin  channels  are 
2R  and  2r,  respectively;  the  net  conductance  is  l/2R+l/2r.  In  (b),  the  resistance 
for  either  channels  are  r  +  R;  the  net  conductance  is  2/(r  + 


1.2.2    Fermi  Surface  Mismatch  Picture 

Although  the  basic  idea  of  the  GMR  can  be  imderstood  by  the  simple  resistor 
model  above,  this  model  does  not  explain  the  origin  of  the  spin  anisotropy,  which 
is  the  difference  divided  by  the  sum  of  the  resistance  in  the  two  spin  channels. 
There  are  two  spin  anisotropics,  the  bulk  spin  anisotropy  and  the  interface  spin 
anisotropy.  The  interface  spin  anisotropy  can  be  understood  with  simple  pictures 


8 

and  is  often  the  dominant  effect  of  the  GMR  unless  the  layers  are  thick.  Therefore, 
in  the  following  we  will  explain  the  GMR  using  only  the  interface  spin  anisotropy. 

The  most  intuitive  way  of  explaining  the  spin  anisotropies  and  the  CPP-GMR 
is  probably  the  Fermi  surface  picture.  We  use  the  Fermi  surfaces  plots  from  our 
"Fermi  Surface  Database"  website  at  http://www.ufl.edu/fermisurface  [SC- 
SI] in  the  following  discussion.  It  provides  both  pictures  and  interactive  Virtual 
Reality  Modeling  Language  (VRML)  models  of  Fermi  surfaces  for  most  elemental 
metals.  We  assume  the  layers  to  be  large  compared  with  the  Fermi  wavelength 
such  that  the  bulk  Fermi  surfaces  is  a  reasonable  approximation.  Again  the  two 
fluid  model  is  used,  assuming  the  spin  diffusion  length  is  large  compared  to  the 
thickness  of  the  layers  [16,  28-29]. 

As  an  example,  we  consider  the  Fe/Cr/Fe  system,  a  three-layered  structure 
formed  by  sandwiching  a  Cr  layer  in  between  two  Fe  layers.  The  Fermi  surfaces 
of  bulk  Fe  spin-up  and  spin-down  channels,  and  non-magnetic  Cr,  are  arranged 
in  Figure  1.6  for  both  (a)  the  parallel  and  (b)  the  anti-parallel  magnetic  config- 
urations. As  shown  in  Figm-e  1.6(a),  the  Cr  Fermi  surface  is  very  different  from 
the  Fe  spin-up  Fermi  siuface,  while  it  is  much  closer  to  the  Fe  spin-down  Fermi 
surface.  If  an  electron  in  the  spin-up  channel  of  the  Fe  layer  on  the  left  needs 
to  go  into  the  Cr  layer,  there  will  be  larger  scattering  at  the  interface  because  a 
wavevector  on  the  spin-up  Fermi  surface  of  bulk  Fe  is  not  on  the  bulk  Cr  Fermi 
surface.  When  this  electron  goes  from  the  Cr  layer  to  the  Fe  layer  on  the  right,  it 
is  scattered  again.  Therefore,  an  electron  goes  though  two  large  scattering  events, 
and  as  a  result,  the  spin-up  channel  has  a  low  conductivity.  On  the  other  hand,  a 
spin-down  electron  from  the  Fe  layer  on  the  left  is  not  scattered  as  much  in  going 
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through  the  two  interfaces.  The  scattering  is  much  smaller  because  the  Cr  and  Fe 
spin-down  Fermi  surfaces  are  similar.  Thus,  the  spin-down  channel  has  a  higher 
conductivity.  The  total  conductivity  for  the  parallel  configuration  is  the  sum  of 
the  two  channels,  which  is  dominated  by  the  spin-down  channel. 


(a) 


(b) 


Spin  Up  Channel 

Cr 

Feup'T 

Spin  Down  Channel 

Fedn 

Cr 

Fedn 

Spin  Up  Channel 

-eup 

Cr 

Fe< 
t 

in 

Spin  Down  Channel 

Fedn 

Cr 

Figure  1.6:  Fermi  Siuface  Mismatch  in  the  Fe/Cr/Fe  multilayer.  The  bulk  Fermi 
for  each  spin  in  the  layers  in  the  (a)  parallel,  and  (b)  anti-parallel,  magnetic  con- 
figuration are  shown.  There  are  only  three  different  Fermi  surfaces,  Fe  majority 
spin,  Fe  minority  spin,  and  Cr.  Since  the  spin  of  an  electron  crossing  the  layers  is 
assumed  to  fixed,  the  electron  behave  as  a  majority  or  minority  electron  depending 
on  whether  its  spin  is  aligned  or  anti-aligned  with  the  local  magnetization  of  the 
layer. 

The  anti-parallel  case  is  shown  in  Figure  1.6(b).  The  spin-up  and  spin-down 
Fermi  smfaces  for  the  Fe  layer  on  the  right  are  exchanged.  When  the  local  magne- 
tization on  the  left  and  the  right  are  anti-parallel  to  each  other,  an  electron  which 
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starts  as  a  spin-up  (majority)  electron  on  the  left  will  be  in  the  minority  population 
on  the  right.  Hence  it  behave  as  a  spin-down  electron  in  the  bulk  of  the  Fe  layer 
on  the  right.  Here,  electrons  in  the  both  spin  channel  go  though  the  same  amoimt 
of  scattering. 

Since  the  Cr/Fe  spin- up  Fermi  surface  mismatch  is  much  larger  than  the  Cr/Fe 
spin-down  Fermi  siurface  mismatch,  the  conductivity  for  the  spin-down  channel 
of  the  parallel  configuration  is  much  larger  than  the  other  channels.  Therefore, 
the  conductivity  for  the  parallel  magnetic  configuration  is  much  larger  than  the 
conductivity  for  the  anti-parallel  magnetic  configuration.  In  this  case,  the  giant 
magnetoresistance  is  positive. 

The  CPP-GMR  in  Co/Cu/Co  can  be  explained  similarly  using  the  Fermi  sur- 
faces in  Figure  1.7.  Here,  the  Cu  Fermi  surface  is  similar  to  the  Co  spin- up  Fermi 
surface,  but  rather  different  from  the  Co  spin-down  Fermi  surface.  Therefore, 
the  conductivity  is  dominated  by  the  spin-up  channel  of  the  parallel  configura- 
tion. Again,  the  conductivity  for  the  parallel  magnetic  configuration  is  significantly 
larger  than  the  conductivity  for  the  anti-parallel  magnetic  configuration.  Hence, 
the  OMR  is  also  positive. 

We  note  that  the  spin  anisotropy  for  the  Fe/Cr  interface  has  a  different  sign 
than  the  spin  anisotropy  for  the  Co/Cu  interface. 

1.2.3    Potential  Step  Picture 

A  related  way  of  understanding  the  spin  anisotropy,  or  the  difference  in  re- 
sistance for  difference  spins,  at  interfaces  is  to  look  at  the  potential  step  across 
the  interfaces.   As  will  be  explained  in  Chapter  2,  the  major  difference  in  the 
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Figure  1.7:  Fermi  Surface  Mismatching  the  Co/Cu/Co  multilayer.  The  bulk  Fermi 
for  each  spin  in  the  layers  in  the  (a)  parallel,  and  (b)  anti-parallel,  magnetic  con- 
figinration  are  shown.  ' 

tight-binding  parameters  in  the  3d  transition  metals  is  the  on-site  energy  of  the 
d-orbitals.  The  system  looks  like  a  simple  transmission-though-potential-barrier 
problem.  We  must  be  careful,  however,  that  the  potential  step  affects  only  the 
d-electrons,  while  the  s-  and  p-electrons  carry  a  significant  fraction  of  the  current. 
The  effect  is  comphcated  by  the  spd-hybridization. 

From  the  tight-binding  parameters,  the  Cu  potential  is  0.82  eV  lower  than  the 
spin-up  Co  potential  and  is  2.3  eV  higher  than  spin-down  Co  potential.  Since 
the  potential  step  in  Cu/Co  spin-up  interface  is  higher,  it  has  higher  resistance. 
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Similarly,  the  Cr  potential  is  1.53  eV  higher  than  the  spin-up  Fe  potential  and  is 
0.23  eV  lower  than  the  spin-down  Fe  potential. 

The  potential  step  may  also  be  used  to  estimate  the  impurity  potential  due  to 
interface  diffusion  of,  say  Cu  into  Co.  ^ 

1.2.4    Effects  of  the  Disorder 

Although  the  simple  pictures  in  the  last  section  can  help  imderstand  the  origin 
of  the  giant  magnetoresistance  effect,  these  pictures  do  not  include  both  the  band 
structure  effect  and  the  disorder  effect.  We  have  explained  that  the  band  structure 
is  often  the  source  of  the  GMR.  On  the  other  hand,  both  the  conductivity  and  the 
size  of  the  GMR  are  determined  by  the  disorder.  If  we  compute  the  conductivity 
with  small  disorder,  we  would  get  a  conductivity  orders  of  magnitudes  larger  than 
what  is  usually  obtained  in  experiments.  Also,  if  the  disorder  is  too  small,  the 
GMR  would  be  too  large  compared  to  experiments.  Therefore,  disorder  must  be 
included  in  the  calculation  in  order  to  obtain  the  conductivity  and  the  GMR. 

As  explained  by  J.  Chen  et  al.  [25],  the  disorder  can  affect  the  the  density 
of  states  (DOS)  significantly.  The  density  of  states  of  Cu  and  Co  with  different 
disorder  parameters  is  shown  in  Figure  1.8.  As  we  can  see,  the  disorder  significantly 
changes  the  density  of  states.  The  effect  of  the  disorder  is  to  spread  out  the  DOS. 
The  net  effect  at  the  Fermi  surface  is  hard  to  predict  without  the  calculations  done 
in  Chapter  2.  In  Cu  and  the  spin-down  channel  of  Co,  the  DOS  at  the  Fermi  level 
is  not  significantly  affected.  However,  for  the  spin-up  channel  of  Co,  which  has  a 
similar  band  structure  to  Cu,  the  DOS  is  doubled  due  to  the  DOS  spread  from 
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the  d-peak  just  below  the  Fermi  level.  We  will  see  in  the  next  chapter  that  this 
significantly  reduces  the  conductivity,  even  though  the  DOS  is  doubled. 


-4.0 


-2.0  0.0 

E-Ep  (eV) 


Figure  1.8:  Density  of  states  (DOS)  of  Cu,  and  Co  spin- up  and  spin-down  channels 
for  different  impmity  density.  When  the  impurity  density  is  small,  the  DOS  remain 
close  to  the  unperturbed  DOS.  However,  when  the  impurity  density  is  moderate, 
the  DOS  can  change  a  lot.  The  conductivity  is  affected  by  the  DOS  near  the  Fermi 
level.  For  Cu  and  Co  spin-down,  the  change  in  DOS  at  the  Fermi  level  are  small. 
However,  for  Co  spin-up  channels,  which  has  a  similar  band  structure  to  Cu,  the 
DOS  at  the  Fermi  level  is  increased  by  a  factor  of  two  [25]. 


To  calculate  the  GMR,  we  have  to  find  the  appropriate  disorder  parameters. 
We  set  the  bulk  impiurity  parameter  so  that  the  bulk  resistivity  is  similar  to  experi- 
ments. A  single  bulk  impurity  parameter  should  be  sufficient  to  give  the  right  bulk 
resistivity  for  all  elements  in  a  sample.  In  Chapter  1,  we  explained  how  to  find 
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the  impurity  potential  due  to  diffusion.  Tlie  density  of  the  impurity  is  harder  to 
set  because  there  is  no  experiment  the  clearly  measure  this  parameter.  However, 
since  only  the  total  scattering  rate  matters,  one  can  set  the  parameters  to  match 
the  experiments.  Since  the  same  impurity  density  has  to  be  used  for  both  parallel 
and  anti-parallel  magnetic  configurations,  the  single  parameter  should  be  able  to 
give  the  correct  resistivity  for  both  configmration.  In  other  words,  there  should  less 
parameter  than  the  niraiber  of  data  points. 

In  Figure  1.9,  the  current-perpendicular-to-plane  resistivity  times  the  layer 
thickness  t  for  the  Co(/Cut  superlattices  is  plotted  against  t.  The  results  for 
the  spin-up  channel,  labeled  as  Cott/Cut,  and  the  spin-down  channel,  labeled 
as  Coit/Cuf,  are  shown  for  calculations  with  bulk  disorder  only,  bulk  plus  spin- 
independent  interface  disorder,  and  bulk  plus  spin-dependent  interface  disorder. 
The  lines  are  Uneax  fits  to  the  computed  results.  The  interface  resistances  can  be 
read  from  the  intercepts.  The  bulk  resistance  is  related  to  the  slope.  Both  the  in- 
terface resistance  and  the  bulk  resistance  are  affected  by  the  types  of  disorder.  It  is 
found  that  both  bulk  disorder  and  spin-dependent  interface  disorder  are  necessary 
to  explain  experiment  results  [25] . 

1.3    Optimization  at  the  Atomic  Level 

Optimizations  in  physical  systems  axe  often  performed  by  the  two  popular 
methods:  simulated  annealing  [32]  which  is  inspired  by  statistical  mechanics,  and 
genetic  algorithms  [33-34],  which  are  inspired  by  evolution.  Simulated  annealing 
has  been  used  in  optimization  problems  such  those  in  circuit  design  [32].  It  is  also 
used  in  problems  such  as  structviral  optimization.  For  example,  it  has  been  used 
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Figure  1.9:  Current-perpendicular-to-plane  resistivity  of  superlattices  Co/Cu.  The 
spin- up  and  spin-down  channels  are  labeled  Cott/Cut  and  CoJ.t/Cut,  respectively. 
The  resistivity  are  calculated  for  three  cases:  with  only  bulk  disorder,  bulk  plus 
spin-independent  interface  disorder,  and  bulk  plus  spin-dependent  interface  disor- 
der. The  lines  are  linear  fits  to  the  computed  results.  The  interface  resistances  can 
be  read  from  the  intercepts  [25]. 

to  find  the  most  stable  structure  of  a  7  Si  atom  cluster  using  the  pseudo-potential 
density-functional  theory  [35].  The  genetic  algorithm  has  also  been  used  in  opti- 
mization of  the  molecular  structure.  For  example,  Deaven  and  Ho  [36]  performed  a 
structural  optimization  by  starting  with  a  population  of  random  molecules  with  60 
caxbon  atoms.  Then  the  more  stable,  or  "fit" ,  molecules  are  cut  into  half  and  glue, 
or  "mate",  with  each  other  to  produce  a  new  population  of  molecules.  Genetic 
manipulations  such  as  crossing  over  are  performed  in  the  process.  As  the  simula- 
tion goes  on,  the  most  stable  structure,  the  ideal  icosahedral  buckyball  structure, 
is  achieved. 

In  Chapter  5,  we  turn  to  a  new  direction  of  research  in  which  we  optimize  the 
atomic  configuration  of  a  nanostructure  to  obtain  the  maximiun  GMR.  Performing 
optimization  of  physical  properties  in  the  atomic  scale  is  an  exciting  new  trend. 
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Pranceschetti  and  Zunger  performed  the  simulated  annealing  optimization  to  find 
the  atomic  configurations  with  the  maximum  bandgap  in  semiconducting  alloys 
and  superlattices  [37] .  To  maximize  the  bandgap  energy  with  simulated  annealing 
(Figure  1.10),  they  started  with  a  random  atomic  configuration  and  performed 
Monte  Carlo  moves  in  the  configuration  space,  which  is  the  space  of  all  possible 
configurations.  A  move  can  be  accepted  with  a  probability  even  when  the  new 
configuration  has  a  smaller  bandgap.  However,  as  the  annealing  goes  on,  the 
probability  is  reduced  until  the  system  finally  rests  in  the  maximum  bandgap 
configuration.  The  maximum  bandgap  configuration  for  a  Alo.25Gao.75As  alloy  is 
shown  in  Figiu-e  1.11.  It  would  be  difficult  to  find  this  configuration  using  trial- 
and-error. 

1-90,  ,  ,  ,  ,  ,  1 


1.74 1  1  1  '  1  1  1 

0    2.000  4,000  6,000  8.000  10,000  12,000 

Monte  Carlo  move 

Figure  1.10:  Optimization  of  the  bandgap  energy  of  Alo.25Gao.75As  alloy  by  sim- 
ulated annealing  [37].  The  bandgap  energy  fluctuates  as  the  Monte  Carlo  rou- 
tine progress.  However,  since  there  is  a  bias  towards  increasing  bandgap  and  the 
amount  of  fluctuation  is  decreasing,  the  system  finally  arrive  a  maximum  bandgap 
configuration. 

Although  atomic  level  optimization  has  been  performed,  the  optimization  of 
physical  properties  other  than  the  energy  is  new.  Here  one  searches  for  an  atomic 
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Figure  1.11:  Maximum  bandgap  configuration  of  Alo.25Gao.75As  [37].  This  compli- 
cated configuration  is  unlikely  to  be  found  by  trial-and-error. 

configuration  that  is,  by  definition,  not  the  most  stable  configuration.  However, 
with  the  non-equilibrium  growth  techniques,  structures  optimized  for  the  desire 
physical  property  may  become  a  reality. 

1.4    Overview  of  the  Rest  of  the  Dissertation 

The  rest  of  the  dissertation  goes  as  follows.  In  Chapter  2,  we  will  develop  a 
method  to  compute  the  conductivity  in  magnetic  nanostructures.  Both  the  band 
structure  and  the  disorder  will  be  included.  In  Chapter  3,  we  will  discuss  how  to 
apply  the  model  efficiently  by  rewriting  the  equations.  In  Chapter  4,  we  apply  the 
method  to  calculate  the  CPP-GMR  and  compare  the  result  with  experiments.  In 
Chapter  5,  we  will  optimize  the  atomic  configuration  of  a  nanostructure  to  obtain 
the  maximum  OMR.  In  the  final  chapter,  we  will  discuss  the  imphcations  of  this 
work  and  future  directions. 


CHAPTER  2 

MODELING  TRANSPORT  IN  NANOSTRUCTURES 


In  this  chapter,  we  discuss  how  to  model  electron  transport  in  magnetic  nanos- 
tructures.  The  common  approaches  include  semi-classical  methods,  simplified  one 
or  two  band  tight-binding  model,  multi-band  tight-binding  model,  and  density- 
functional  calculations.  Each  model  contains  an  increasing  level  of  details  and  has 
its  advantages  and  disadvantages. 

Semi-classical  methods  [16]  can  be  used  for  studying  large  systems  with  complex 
geometry  that  other  methods  can  not  solve  efficiently.  The  conductivity  is  often 
calculated  with  the  Boltzmann  equation.  Bulk  and  interface  scattering  effects 
can  be  included  as  parameters  in  the  Boltzmann  equation,  without  the  need  for 
microscopic  details.  This  approach  valid  when  the  length  of  smallest  feature  is 
much  larger  than  the  Fermi  wavelengths,  but  it  is  often  possible  to  extend  it  to 
smaller  systems  for  qualitative  study.  However,  when  the  feature  length  scale  is 
only  a  few  nanometers,  the  semi-classical  method  may  become  problematic. 

The  simplified  one  or  two  band  tight-binding  model  is  often  used  to  give  qual- 
itative understanding  of  the  finite  size  effects  and  some  of  the  band  structiure 
effects.  Usually,  the  parameters  in  this  model  axe  obtained  by  reasonable  esti- 
mations. The  conductivity  is  usually  obtained  by  the  Kubo  formula  [38]  or  the 
Landauer  formula  [39-40].  This  model  is  easy  to  solve,  and  sometimes  it  is  even 
possible  to  obtain  analytical  solutions.  In  addition,  the  simplified  band  structure  is 
helpful  in  understanding  some  effects  that  may  be  less  obvious  in  more  complicated 
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models.  However,  in  real  materials,  such  as  iron  or  nickel,  the  band  structxure  is 
too  complicated  for  the  results  of  this  simple  model  to  compare  with  experiments 
quantitatively. 

The  multi-band  tight-binding  model  [22-26]  is  used  for  quantitative  compar- 
ison with  experiments.  Parameters  are  often  fit  from  density-functional  theory 
calculations  and  hence  it  captures  the  band  structure.  Details  of  atomic  varia- 
tion, impurity  scattering  effects  [25-26]  can  be  modeled  accurately.  Although  this 
method  is  not  as  accurate  as  a  full  density-function  theory  calculation,  it  has  the 
advantage  of  being  much  faster.  It  is  also  slower  than  the  simple  tight-binding 
model;  however,  as  we  will  show  in  detail  in  Chapter  3,  it  is  fast  enough  to  per- 
form large  scale  studies  in  inexpensive  modern  computer  platforms  such  as  a  small 
Pentium  III  cluster. 

The  density- functional  theory  is  the  most  accm-ate  method  [27,  41-42],  and 
unlike  the  above  mentioned  methods,  less  insight,  or  educated  guess,  is  required 
in  the  calculations.  The  density-functional  calculations  can  also  provide  informa- 
tion for  the  less  demanding  tight-binding  calculations,  and  is  useful  for  vaUdat- 
ing  the  results  obtained  from  other  studies.  However,  it  demands  much  greater 
computer  resources  than  tight-binding  methods.  Recently,  large  scale  density- 
fimctional  studies,  such  as  the  complex  magnetic  properties  [43],  are  performed 
with  the  massively  parallel  processing  computers. 

In  the  following  sections,  we  first  discuss  the  multi-band  tight-binding  Hamil- 
tonian  and  the  impiurity  averaged  Green's  function.  Then  we  will  explain  the 
calculation  of  the  conductivity  using  the  Kubo  formula  and  the  vertex  correction 
needed  for  current  conservation.   In  the  last  section,  we  discuss  how  to  obtain 
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tight-binding  and  impurity  parameters  corresponding  to  experimental  situation 
and  compare  some  of  the  theoretical  results  with  experiments. 

2.1    Modeling  of  the  Band  Structure 

In  this  section,  we  discuss  the  modeling  of  the  band  structure  using  a  tight- 
binding  Hamiltonian.  Then  we  explain  how  to  model  impiurities  in  samples  and 
how  to  obtain  the  impurity  averaged  Green's  function,  which  contains  the  band 
structure  information. 


2.1.1    Multi-band  Tight-binding  Hamiltonian  . 

The  general  form  of  a  tight-binding  Hamiltonian  is  given  by 

H=  H{r,r';j,j';a,a')clj,Cr'ry,  ..  (2.1) 

where  clj^  and  Crj>  are  creation  and  annihilation  operators  of  an  orbital  j  with 
spin  a  at  site  r.  Depending  on  the  crystal  structure  and  the  details  required  in  the 
model,  the  hoping  terms,  H{r,T';  jj';a,a'),  are  set  to  be  zero  when  r  and  r'  are 
third  neighbors  or  fiu-ther.  However,  in  bcc  crystals,  the  third  neighbor  is  essential. 
In  modeling  the  transition  metals,  it  is  often  important  to  keep  all  the  orbitals  in 
the  outermost  shell:  s,p,d. 

The  coupling  between  the  majority  and  minority  spins  H{r,r';j,j';^,i)  con- 
tributes to  spin-flip  scattering.  For  example,  to  include  the  spin-orbit  coupling,  a 
term  •  S  is  added  to  the  Hamiltonian  [44].  For  transports  in  transition  metals 
samples  we  consider,  the  spin-diffusion  length  is  much  larger  than  the  layer  thick- 
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ness  [28-29],  therefore  spin-flip  is  neglected  in  most  calculations.  In  this  case,  the 
two  spin  species  axe  assumed  to  be  independent  and  with  decoupled  Hamiltonians, 
i/(t,  T)  ^■nd  H{l,i),  and  can  be  solved  separately. 

Unless  otherwise  specified,  the  formulation  is  this  work  can  be  applied  to  the 
full  Hamiltonian  as  well  as  the  Hamiltonian  for  individual  spins. 

2.1.2    Modeling  Impurities 

Scattering  sources  exist  in  all  samples  including  in  experiments  of  heterostruc- 
tiu-es.  If  we  do  not  include  these  effects  in  our  calculations,  the  resistivity  computed 
would  be  at  least  one  or  two  magnitudes  less  than  those  observed  in  experiments. 
There  are  three  major  source  of  scattering:  disorder  in  the  bulk  of  the  sample, 
impurity  atoms  usually  due  to  interface  diffusion,  and  phonon  scattering. 

Figure  2.1  shows  a  clean  multilayer  lattice  and  the  same  lattice  with  interface 
diffusion,  where  10%  of  each  monolayer  diffuses  into  each  of  the  two  monolayers 
above  and  below.  One  can  think  of  an  impiurity  atom  in  a  host  as  a  new  het- 
erostructure.  However,  as  we  will  explain  in  the  next  subsection,  the  advantage  of 
treating  it  as  an  impurity  is  that  simplification  is  often  possible.  For  example,  in 
an  Co/Cu  interface,  we  would  treat  the  Co  atoms  diffused  into  the  Cu  layers  as  Co 
impurity  in  Cu  host,  and  the  Cu  atoms  diffused  into  the  Co  layers  as  Cu  impurity 
in  Co  host.  The  impurity  potential  of  an  impurity  B  in  an  host  A  is  defined  as 


A«(r,r')   =  HB{r,v')-HA{r,r') 


(2.2) 
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where  Hb{t,  r')  and  Ha{t,  r')  are  the  hoping  between  r  and  r'  for  the  perturbed  and 
unperturbed  system.  The  most  important  contribution  to  Aw  is  usually  the  on-site 
term.  As  a  perturbation,  we  found  it  sufficient  to  model  this  term  as  diagonal  in 
all  indices.  In  other  words, 

Au{r,T';j,j';a,a')  =  Au{r,T-j,j;a,a)5rr'Sjj'S^a'.  (2.3) 


Figure  2.1:  Illustration  of  a  multilayer  without  and  with  interface  diffusion.  With- 
out the  interface  diffusion  (a),  the  multilayer  has  translational  symmetry  in  the 
in-plane  directions.  In  (b),  10%  of  each  monolayer  diffuse  into  each  of  the  two 
monolayers  above  and  below.  By  adding  interface  diffusion,  the  translational  sym- 
metry is  broken. 

Two  types  of  disorders  are  included  in  our  model:  the  substitutional  dis- 
order which  is  often  caused  by  interface  diffusion  and  the  structural  disorder 
which  are  randomly  located  in  the  sample.  To  model  the  substitutional  dis- 
order in  the  3d  transition  metals,  we  set  the  substitutional  disorder  parameter 
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Awsub(r,  r;j,j;  a,a)  =  0  for  j  —  s,Px,Py,Pz  because  the  most  important  difference 
between  the  tight-binding  parameters  in  different  elements  is  the  on-site  poten- 
tial of  the  d-orbitals.  The  difference  among  the  d-orbital  parameters  are  small; 
therefore,  it  is  often  a  good  approximation  to  set  Attsub(r>  r;  ifi;  cr)  =  AMd(r,  cr) 
for  j  =  dxy,dyz,dzx,di,d2,  such  that  two  parameters,  AMd(r,  j)  and  Awd(r,  J,), 
parameterize  the  impurity  at  site  r. 

To  model  structural  disorders,  we  use  impurity  potential  of  the  form 

Au(r,  r';  j,  f;  a,  a')  =  AustSrr'SjfS^a' ,  (2.4) 

where  Aitgt  is  a  single  parameter  that  characterizes  the  bulk  of  the  sample.  In 
this  work,  we  do  not  model  the  phonon  scattering.  However,  we  can  increase  the 
structural  disorder  parameter  for  the  bulk  to  get  an  idea  of  the  effect  of  increasing 
temperature. 

2.1.3    Impurity  Averaged  Green's  Functions 

Once  the  Hamiltonian  and  the  impurity  parameters  are  given,  we  are  ready 
to  obtain  the  band  structure  information  by  either  solving  the  eigenvalue  problem 
or  the  Green's  fvmction.  Here  we  explain  the  impurity  average  Green's  function 
approach  which  is  especially  useful  in  treating  systems  with  impurity  scattering. 
The  impurity  average  Green's  function  is  the  ensemble  average  of  the  Green's 
function  for  systems  with  different  microscopic  impmrity  configiurations  in  the  same 
macroscopic  impurity  distribution.  In  many  systems,  geometric  symmetries  are 
broken  by  the  impiu-ity  atoms.  For  example,  in  the  multilayer  structiire  in  Figmre 


24 

2.1(b),  the  in- plane  translational  symmetries  in  the  unperturbed  system  (Figure 
2.1(a))  are  broken  by  the  impurities  due  to  interlayer  diffusion.  However,  the 
impurity  probabihty  distributions  in  each  site  within  the  same  layer  are  the  same. 
Therefore,  after  averaging  over  the  impurity  ensemble,  the  system  recovers  the 
translational  symmetry.  Also,  the  average  of  the  Green's  function  is  performed 
analytically.  As  a  result,  the  impurity  averaged  Green's  function  approach  can  be 
used  to  study  much  larger  systems  than  techniques  that  requires  averaging  over 
the  impurity  configuration  numerically. 

We  first  consider  the  system  unperturbed  by  impurities.  The  unperturbed,  or 
bare,  retarded/advanced  Green's  functions  is  given  by 

{uj-  H±ir])Go^'^  =  1,  ■       "  .  •  (2.5) 

where  the  scalar  u  is  the  energy,  t]  is  an  infinitesimal  small  positive  number,  and 
1  is  the  identity  operator.  The  Green's  function  and  the  quantities  obtained  from 
the  Green's  function  are  a;-dependent,  in  the  following,  the  u)  label  is  suppressed 
for  clarity.  From  the  definition,  the  bare  Green's  functions  satisfy  : 

{G^^y  =  Gi^^.  ,  (2.6) 

Since  the  tight-binding  parameters  we  use  are  real,  the  Hamiltonian  is  real^  in  the 
real  space  representation.  In  other  words,  H'^  =  W  —  H,  which  leads  to 


This  is  not  true  in  the  Fourier  space 


(2.7) 
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When  impurity  atoms  are  added  to  the  system,  one  may  try  to  compute  the 
Green's  function  for  system  with  impurities  at  their  exact  location.  However,  in 
most  experiments,  the  exact  location  of  impurity  atoms  are  unknown,  and  even  the 
probability  distribution  of  the  impurities  are  obtained  by  educated  guesses.  More- 
over, the  effect  of  the  randomly  distributed  impurities  are  averaged  out.  Therefore, 
one  has  to  compute  the  impurity  averaged  Green's  fimction^,  defined  as 

{u-H-  S«/^)G^/^  =  1.  (2.8) 

where  S^/"^  is  the  self-energy.  For  the  on-site,  diagonal  impurity  potential  de- 
scribed in  the  last  section,  the  self-energy  is  given  by 

S«/^(ri,r2)   ^  ^(ri)Jr,r,G'«/^(ri,r2)  (2.9) 
where  the  impurity  parameter 

^(ri)   =   {Au\r^))  (2.10) 

is  the  ensemble  average  of  the  square  of  the  impurity  potential.  For  convenience, 
we  denote  Equation  2.9  as 

Y^R/A  ^  j;^^.  (2.11) 

^For  convenience,  we  use  G^/-^  instead  of  {G^/^)  as  the  symbol  for  the  impurity  averaged 
Green's  function. 
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In  this  formula,  we  have  neglected  the  cross-diagrams,  which  corresponds  to  fourth 
order  terms  in  the  impurity  perturbation  [45-46].  Since  the  impurity  average 
Green's  function  and  the  self-energy  depend  on  each  other,  we  have  to  solve  Equa- 
tions 2.8  and  2.9  self-consistently.  Other  types  of  scattering  effects,  such  as  those 
by  the  bulk  disorder,  can  be  modeled  within  the  same  framework. 
We  can  also  write  the  Green's  functions  in  Dyson's  equation  as 

G^M  =  G?/^  +  gJ/'^E^/^G^/^.  (2.12) 

The  following  properties  follows  from  Equation  2.6: 

(G«/^)t  =  G^/^ 

(E«M)t  =  s^/fi.  (2.13) 

In  addition,  in  real  space  representation,  ,  . 

(S«/^)^  =  S«/^.  ';  (2.14) 

2.2    Conductivity  Calculation 

In  this  section,  we  explain  how  to  compute  the  conductivity  using  the  Kubo 
formula  which  is  the  theory  of  linear  response.  We  first  review  different  current 
related  operators  in  the  tight-binding  model,  which  will  be  then  be  used  in  the 
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Kubo  formula.  We  will  also  discuss  the  vertex  correction  needed  to  give  a  current 
conserving  conductivity. 

2.2.1    Formal  Aspects  of  Current  Operators 

In  this  and  the  next  subsections,  we  discuss  discuss  four  operators  related  to  the 
current:  the  total  current,  the  divergence  of  current  density,  the  current  through  a 
siuface,  and  the  current  density  operators.  These  operators  will  be  useful  for  com- 
puting the  conductivity  and  checking  the  ciurrent  conservation.  In  this  subsection, 
we  will  discuss  these  operators  formally.  In  the  next  subsection,  we  express  these 
operators  in  matrix  form  explicitly  for  an  one  dimensional  Hamiltonian. 

Classically,  the  current  density  is  defined  as  the  charge  density  times  the  veloc- 
ity. The  total  current,  the  divergence  of  current  density,  and  the  current  through  a 
surface  can  be  obtained  from  integrating  or  differentiating  the  current  density.  In 
quantum  mechanics,  the  current  operator  is  defined  using  the  Hamiltonian.  How- 
ever, in  the  tight-binding  model,  the  current  density  is  more  difficult  to  define 
because  of  the  discrete  nature  of  the  basis. 

The  current  density  represented  in  a  discrete  basis,  it  is  defined  from  the  current 
from  atomic  sites  r  to  r'  divided  by  the  cross-sectional  area.  In  other  words,  there  is 
only  off  diagonal  element  in  the  matrix  representation.  Unfike  the  classical  current 
density,  the  cvurent  density  operator  not  a  localized  in  space.  On  the  other  hand, 
the  total  current,  the  divergence  of  current  density,  and  the  current  operators  are 
easily  found  from  the  charge  density  operators  though  the  current  conservation 
equations  in  Table  2.1.  We  will  use  the  current  operator  to  define  the  current 
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name 

operator 

relation  to  Jr 

conservation  equation 

total  current 

J 

divergence 

Vr  - Jr 

^a=l  aa 

current 

/ 

an 

^  —  ~'§i  S(r-ro)  n<0  ^rPr 

Table  2.1:  Summary  of  current  density  related  operators  in  the  tight-binding 
model.  The  names,  symbols,  relations  to  the  current  density  operator  Jr,  and 
conservation  equations  used  are  shown  in  the  table.  In  the  tight-binding  model, 
we  define  the  current  operators  using  the  conservation  equations,  and  use  them  to 
define  the  current  density  operator.  In  the  table,  Vr  is  the  volume  per  atom,  pr  is 
the  charge  density,  P  is  the  polarization  operator,  is  the  atomic  spacing  in  the 
Qf-th  direction. 

density  operator,  and  use  the  operator  for  the  total  current  and  the  divergence  of 
the  cmrrent  density  to  check  om-  definition. 

The  relations  between  the  operators  discussed  in  this  section  are  similar  to 
those  of  the  classical  quantities;  however,  some  details  must  be  treated  carefully. 
Table  2.1  summaries  the  current  operators  and  how  can  they  be  obtained  thru 
current  conservation  laws. 

The  total  cmrent  operator  can  be  obtained  [47]  from  the  polarization  operator, 
which  is  defined  as 


P  =  e  ^  rc|:cr. 


(2.15) 
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The  total  current  operator  is  given  by 

J   =  ^ 

dt 

=  i[H,P] 

=  -ie^{T-T')H{r,r')clcr'.  (2.16) 

r,r' 

Here,  the  Planck's  constant  divided  by  27r,  is  set  to  equal  to  one.  Since 
classically  the  total  current  is  the  volume  integral  of  the  current  density,  one 
may  be  tempted  to  define  the  current  density  operator  Jr  as  proportional  to 
X^r'(l/K')(r  ~  r')[/f(r, r')4cr'  —  h.c.].  However,  this  definition  violates  current 
conservation  because  the  divergence  of  this  expression  does  not  agree  with  the 
-dpr/dt. 

We  find  the  operator  Vr  •  Jr  using  the  continuity  equation 

K.J,  =  -^.  (2.17) 
In  tight-binding  model,  the  density  is  given  by 

Pr  =  ec\CrlVr,  (2.18) 

where  Vj.  is  the  volume  occupied  by  the  atom  at  site  r.  The  Schrodinger  equation 
implies 


Vr-Jr  = 


—ie 


Vr 


[Hclcr  -  clcj-H] 


(2.19) 
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—le 


Y,[H{r\r)clcr-clcr'H{r,r')]. 


(2.20) 


It  is  useful  to  note  that  in  matrix  representation  the  only  nonzero  element  of  4cr' 
is  1  at  (r, r').  In  other  words,  (4cr')(ri, =  (^rri(^r'r2-  Therefore,  we  can  write 
the  element  of  the  divergence  of  the  current  explicitly  as 


(Vr  •  Jr)  (ri,r2)  -  i-ie/V,)[H{r,,r)Srr,-5r,rH{r,r2)] 


(2.21) 


The  operator  for  a  current  through  a  plane  containing  point  Tq  and  normal  to 
a  unit  vector  n  can  be  found  by 


(r-ro)  n<0 


=  —te 


5]F(ri,r2)4^Cr2  ,  4cr 

'■ir2  (r-ro)-n<0 

E  E[^(ri>r)4iCr-/i.c.] 

(r-ro)  n<0  ri 


-^^1     E         E     +     E         E      ]  {H{vur)cl^cr  -  h.c] 

(r-ro)-n<0  (ri-ro)-n>0     (r-ro)-n<0  (ri-ro)-n<0/ 

=   -i^     E         E     [H{r,,r)cl^cr-h.c.].  (2.22) 

(r-ro)-n<0  (ri-ro)  n>0 

We  have  used  the  fact  that  the  second  term  in  the  fourth  line  is  zero  due  to  the 
symmetry  between  r  and  ri.  Physically,  the  second  term  is  zero  because  it  contains 
only  the  hoping  term  that  does  not  leave  one  side  of  the  plane.  The  plane  current 
operator  obtained  at  the  end  contains  only  the  hopping  terms  which  bring  current 
from  one  side  of  the  plane  to  another. 
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Unlike  the  operators  shown  in  the  first  column  of  Table  2.1,  there  is  no  simple 
equation  that  leads  to  the  current  density  operator.  In  the  continuous  limit,  the 
current  density  at  r  is  given  by  the  current  though  a  infinitesimal  siurface  divided 
by  the  area.  However,  in  the  tight-binding  model,  we  have  to  discretize  the  length 
scale  to  the  atomic  scale.  As  an  example,  in  one  dimension  the  current  density 
operator  in  second  quantized  form  is  given  by 

=   -ia/K  [H{xi,X2)cl^c^^-  h.c.].    "  (2.23) 

Il>I>I2 

One  can  check  that  the  divergence,  defined  in  one  dimension  as  {Jx+i  —  Jx)/0', 
agrees  with  Equation  2.27,  and  J^Vx  equals  the  total  cm-rent  J.  This  definition 
can  be  extended  to  higher  dimensional  lattices.  For  example,  in  a  rectangular 
lattice  with  only  nearest  neighbors,  the  two  component  of  the  current  density  is 
given  by 

=   -iai/Vr  [H{ri,T2)cl^Cr^-  h.c], 

xi>x>x2,yi=y2=y 

Jr'  =   -ia^/Vr        Yl        [H{rur2)clcr,-h.c.],  (2.24) 

yi>y>y2,xi=x2=x 

where  ai  and  02  are  the  lattice  spacings  in  the  x-  and  ^-directions  respectively. 

2.2.2    Current  Operators  in  One  Dimension 

In  this  subsection,  we  consider  an  one  dimensional  Hamiltonian  and  express 
the  corresponding  current  operators  in  matrix  form  exphcitly.  We  can  better  un- 
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derstand  the  properties  of  the  current  operators  by  studying  this  one  dimensional 
example. 

Consider  a  one  dimensional  Hamiltonian,  H,  with  up  to  second  neighbor  hoping 
terms^.  In  matrix  form,  the  Hamiltonian  can  be  written  as  The  Hamiltonian  in 
matrix  from  is  given  by 


/ 


H  = 


\ 


H{b,3)  H{5A)  H{5,5)  H{5,6)  H{5J) 

if (6, 4)   H{6,5)   H{6,6)   H{6,7)  H{6,8) 

H{7,b)  H{7,6)   H{7,7)   H{7,8)  H{7,9) 


\ 


J 

(2.25) 


From  Equation  2.16,  the  total  current  operator  J  can  be  found  by  —iea  times  the 
matrix 


\ 


2H{5,3)    H{5A)        0  -i/(5,6) 
2//(6,4)    H{6,5)  0 

2i/(7,5)  H{7,6) 


-2H{5,7) 

-H{6J)  -2H{6,8) 

0         -H{7,S)  -2H{7,9) 


V 


where  a  is  the  lattice  constant. 


^The  model  with  only  first  nearest  neighbor  is  too  simple  to  illustrate  some  of  the  properties. 
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The  matrix  form  of  the  current  density  operator  at,  for  example,  x  =  6,  can  be 
obtained  from  Equation  2.23  as 


—tea 


\ 


-i/(5,6)  -F(5,7) 


i/(6,4)  /f(6,5) 
^(7,5) 


(2.26) 


where  a  is  the  lattice  constant,  and  V  is  the  volume  occupied  by  a  site,  using  this 
definition,  one  can  see  that  the  total  current  operator  satisfies  J  =  Ex  JxV. 

We  also  note  that  in  one  dimension,  the  current  operator  is  equals  the  current 
density  times  the  area  V/a.  In  other  words,  the  current  operator  at  x  =  6  contain 
only  the  hoping  terms  which  bring  the  current  from  one  side  of  the  plane  x  =  6" 
to  the  other  side. 

By  Equation  2.26,  the  divergence  of  the  current  density  operator  at  x  =  6  is 
given  by 


dx 


x=6 

a 
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-le 
T 


-H{6,4)   -H{6,b)  0 

H{7, 7) 


-H{6,7)  -H{6,8) 


(2.27) 


This  equation  can  be  rewritten  in  the  following  form: 


dx 


-le 


x=6 


V 


Hi4,6) 

0    0  H{6,6)    0  0 
H{7,e) 
H{8, 6) 
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/ 


0 


0 


//(6,4)  Hi6,5)   ff(6,6)   H{6,7)  H{6,8) 


0 


0 


(2.28) 


which  the  same  as  Equation  2.20  in  one  dimension.  If  we  sandwich  this  expression 
by  two  Green's  functions  or  wavevectors,  the  first  and  the  second  terms  of  the 
expression  can  operate  on  the  left  and  right  respectively  like  the  Hamiltonian.  For 
example,  the  second  matrix  operated  on  the  retarded  Green's  function  on  the  right 
can  be  evaluated  using  Equation  2.8: 


This  property  of  divergence  of  the  current  density  is  useful  later  in  the  disciission 
of  current  conservation. 


5{x,  x). 


(2.29) 
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2.2.3    Kubo  Formula  .  . 

In  the  linear  response  theory,  the  non-local  conductivity  tensor  is  defined  as 
the  response  function  in  the  equation  ,      ;  ' 

Jr  =  5:a„'Er'K',  (2.30) 
r' 

where  Jr  is  the  current  density  response  at  r,  Er'  is  the  applied  electric  field  at  r', 
and  Vf'  is  the  volume  per  atom.  The  impurity  averaged  conductivity  tensor  can 
be  computed  by  the  Kubo  formula  : 


1  rr, 

47r 


(2.31) 


In  these  equations,  the  ciurent  density  plus  vertex  correction  T^^  is  given  by  the 
Dyson's  equations 

r;*^  =  Jr  +  K^^^  (2.32) 


=  3r  +  vG^r^^G^  (2.33) 


=  Jr  +  vG^JrG^  +  vG^vG^3rG^G^  +  ....  (2.34) 

The  other  vertex  corrections  K^^,  K^^,  and  K^^  are  obtained  by  replacing  the 
corresponding  indices  in  the  equations  above.  The  Kubo  formula  is  represented 
in  the  Feynman  diagram  is  Figure  2.2.  Since,  K^'^  is  diagonal  and  F^^  is  not,  in 
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practice,  it  is  more  convenient  to  iterate  on  K^^  instead  of  F^^: 


(2.35) 


vG^'irG'^  +  vG^'K^^G^ . 


(2.36) 


We  also  define 


r' 


Air     y  ^  ^  ^  \ 


(2.37) 


For  a  constant  electric  field,  Jr  =  arE.  Therefore  Oj.  is  often  compared  to  experi- 
ments. It  can  be  shown  that  this  conductivity  is  still  correct  even  in  non-uniform 
electric  field  because  the  result  depends  only  on  the  potential  drop  across  the 
sample. 

As  seen  in  the  bubble  diagrams  in  Figure  2.2,  the  conductivity  can  be  found 
by  first  dressing  up  the  vertex  in  either  current  operator  and  then  multiplying  it 
by  the  other  current  operator,  followed  by  taking  the  trace.  Analogous  to  J,  we 


define  T^^  =  ErTf^K  and        =  Zr^^^Vr,  which  gives 


=  3  +  K 


(2.38) 


=  J  +  i;G^r^^G^ 


(2.39) 


J  +  ?;G^  JG^  +  vG^  vG^JG^  G^  + 


(2.40) 
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Similarly,  the  iteration  equation  for  K^^  is 


(2.41) 
(2.42) 


Therefore,  we  also  have 


2tRAqAj^qR  _  Y^^G^J^G^  -  T^^G^J.G^ 


(2.43) 


As  we  will  discuss  in  the  next  chapter,  the  two  sets  of  equations.  Equations  2.32- 
2.37  and  2.38-2.43,  have  different  computational  costs  in  different  situations.  In 
Chapter  3,  we  will  explain  which  set  of  equations  should  be  used  in  cvurent- 
perpendicular-to-plane  and  current-in-plane  calculations. 


-I- 


Figure  2.2:  Feynman  diagrams  for  the  average  conductivity,  (a)  The  open  bubble, 
which  corresponds  to  the  conductivity  without  the  vertex  correction,  (b)  The 
conductivity  contribution  from  the  vertex  correction. 
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2.2.4    Current  Conservation 

For  a  system  without  impiirity  average,  the  conductivity  is  given  by  the  the 
Kubo  formula  without  the  vertex  correction,  or  the  open  bubble  (Figure  2.2(a)). 
However,  when  averaged  over  the  impurity  configurations,  the  open  bubble  does 
not  satisfies  current  conservation.  This  is  because  using  the  impurity  averaged 
Green's  function  to  compute  the  conductivity  in  the  open  bubble,  one  neglects  the 
correlation  terms  such  as  those  in  Figure  2.2(b).  These  terms  are  required  to  give 
a  conserving  current.  In  the  following,  we  show  the  terms  in  the  vertex  correction 
shown  in  Equation  2.34  will  give  a  conserving  current.  In  other  words,  we  will 
show  that  the  conductivity  obtained  from  Equation  2.37  satisfies 

Vr-Jr     =     Vr-(7rE-0.  (2.44) 

Since  the  electric  field  is  arbitrary,  we  need  to  show 

Vr-(Tr    =    0.  (2.45) 

Since  the  current  Jr'  does  not  depend  on  r,  we  can  first  find  the  divergence  of 
the  vertex  G^FrG^,  and  then  multiply  the  divergence  with  the  Jr',  followed  by 
taking  the  trace.  The  divergence  of  the  zeroth  order  term  of  the  AR  vertex  in 
Equation  2.34  is 

Vr-(G^JrG'«)     =    G^{Vr-Jr)G''.  (2.46) 


In  the  one  dimensional  case,  this  is  simply 


a  ^  '  a 
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(2.47) 


By  current  conservation  (Equation  2.19),  we  have 


(Vr  •  Jr) 

=   {-ie/Vr)[{G^H)  4crG^  -  G^lcr  [HG"") 
=   (-ie/Vr)[(-l  +  G^uj  -  G^S'^)4crG^  -  G\lcr{-1  +  uG""  -  S^G^)] 
=   {-ielVr)[G^clcr  -  4crG^  -  G-^E^cJcrG^  +  G^4crE^G^].  (2.48) 


Here,  we  have  used  the  Schrodinger  equation  (2.8)  to  evaluate  HG^/^  as 


(2.49) 


For  the  divergence  of  the  first  order  term  of  the  AR  vertex  (vertex  with  one 
rung)  in  Equation  2.34,  we  have 


Vr-(G^i;(G^JrG^)G^) 
G^'t;(G^(Vr  •  Jr)G^)G^ 


G^  vG^clcr  G«  -  G^  vclcrG''  G^ 


-G^  vG^i:\lcrG''  G^  +  G^  vG^lcr^'^G''  G^ 


(2.50) 


■  ■'5 
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The  first  term  of  the  above  equation  can  be  rewritten  using  the  following  equation: 

(t;G^4cr)(ri,r2)   =  v{ri)§r,,v2{G^4cr.){ri,r2) 

=  'y(ri)5ri,r2G^(ri,r)(^r,r2 

=  S^(ri,r)5r,r2 

=  (S^4cr)(ri,r2).  (2.51) 

We  can  see  the  first  term  of  Equation  2.50  cancels  the  third  term  of  Equation  2.48. 
Similarly,  the  second  term  of  Equation  2.50  cancels  the  last  term  of  Equation  2.48. 
The  same  argument  can  be  applied  for  higher  order  terms  to  show  the  last  two 
terms  of  the  divergence  of  the  (n  +  l)-th  order  vertex  cancel  the  last  two  terms 
of  the  divergence  of  the  n-th  order  vertex.  As  a  result,  the  only  terms  left  of  the 
divergence  of  the  vertex  are  the  first  two  terms  of  Equation  2.48.  The  divergence 
of  the  conductivity  in  Equation  2.37  is  proportional  to  the  sum  of  the  following 
terms 

2IV(JG^Vr  •  (r;^^)G^)   =   {-2ie/Vr)[Kf  -  4], 
-Tr(JG«Vr  •  (rr^«)G^)   =  {ze/V,)K-4], 
-Tv{3G^Vr-{T^^)G^)   =   {ie/Vr)[Kf  -  4],  (2.52) 

where 


Kf/^  =  TV(JctcrG^/^)  =  Tr(4crG^/^J), 
kJ/^   ^  Tr(JG«/^4cr), 


(2.53) 


In  the  real  space  representation,       =  H,  which  impUes 

(4crG«/^J)^  =  -JG«/^4cr.  (2.54) 

Since  transpose  of  a  matrix  does  not  affect  its  trace,  we  have 

/.f/^  (2.55) 

Unhke  the  continuous  case  [46],  the  divergences  of  the  individual  dressed  bubbles 
in  the  tight-binding  model  are  non-zero.  However,  they  cancel  each  other  so  that 
the  divergence  of  Gr  obtained  from  the  Kubo  formula  is  zero.  Therefore  the  vertex 
corrected  conductivity  satisfies  current  conservation. 

In  this  work,  we  only  use  on-site  scattering  potentials,  the  self-energy  E^/'^  is 
also  on-site.  Therefore, 

E^/^4cr  =  4crS^/^.  (2.56) 

This  means  that  in  Equation  2.48  the  last  two  term  for  divergence  of  the  AA 
and  RR  vertex  cancels  each  other  out,  and  no  vertex  correction  is  needed  for  the 
AA  and  RR  vertex.  The  Kubo  formula  for  on-site  scattering  potentials  can  be 
simplified  as 

ar  =   ^Tr[2JG'^(Jr  +  K,^^)G'^-JG«JrG'^- JG^J.G^],  (2.57) 
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When  the  conductivity  is  calculated  in  nanostructures,  for  example  in  the 
current-perpendicular-to-plane  (CPP)  magnetic  multi-layers,  we  find  the  contribu- 
tion from  the  vertex  corrections  can  be  as  large  as  20%  of  the  total  conductivity. 
Therefore,  the  vertex  correction  should  not  be  neglected. 

2.3    Realistic  Parameters 

In  the  earlier  part  of  this  chapter,  we  explained  how  to  obtain  the  conductivity 
once  the  Hamiltonian  and  the  impurity  potential  are  given.  In  this  section,  we 
explain  how  to  obtain  the  tight-binding  parameters  and  the  impurity  potential, 
which  are  essential  to  modeling  the  electron  transport  accurately.  We  will  also 
compare  some  theoretical  results  with  experiments. 

2.3.1    Tight-Binding  Parameters 

The  tight-binding  parameters  of  an  element  is  obtained  from  fits  to  the  bulk 
band  structure  of  density-functional  calculations  [48].  To  do  this,  we  compute  the 
eigenvalues  for  a  set  of  k-points  in  the  irreducible  zone  with  the  density-functional 
theory,  and  also  with  the  tight-binding  model.  The  tight-binding  parameter  are 
then  adjusted  to  minimize  the  error 

EENf(k)-^f^(k)P,  (2.58) 
k  j 

where  ujj^{k),  and  u;?*^'^(k)  are  respectively  the  tight-binding  theory  eigenvalue 
and  density-functional  theory  eigenvalue  of  the  j-th  orbital  of  a  k- vector.  As  shown 
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by  Papaconstopoulos  [48],  tight-binding  parameters  obtained  in  this  way  give  the 
correct  band  structure  for  bulk  systems. 

For  heterostructures,  we  do  not  obtain  the  tight-binding  parameters  by  fitting 
to  the  band  structure  of  the  heterostructure  because  the  density-functional  theory 
for  complicated  systems  can  be  expensive  computationally.  Since  our  goal  is  to 
model  the  system  efficiently  but  still  give  reliable  predictions,  we  use  a  less  costly 
approach.  The  band  structures  for  3d  transition  metals  are  very  similar,  except  for 
the  offset  of  the  c?-bands,  which  are  different  in  each  metal.  To  model  3c?  transition 
metal  heterostructures,  we  use  the  bulk  parameters  in  the  bulk  and  the  interface. 
The  on-site  energy  of  each  metal  is  shifted  so  that  the  Fermi-level  level  is  zero. 
In  this  case,  the  Fermi  level  for  different  metal  aligns.  The  hoping  parameters 
between  the  atoms  of  two  difference  element  are  set  as  the  average  of  their  bulk 
values. 

2.3,2    Impurity  Potential  for  Substitutional  Disorder 

As  explained  in  Section  2.1.2,  the  substitutional  disorder  is  caused  by  impm-ity 
atoms  usually  from  interface  diffusion.  The  impmrity  potential  for  the  atom  can 
be  obtained  from  density-functional  fits.  To  the  first  order,  one  may  use  the 
difference  in  the  bulk  on-site  parameters  between  the  impurity  and  the  host  as 
the  substitutional  disorder  potential,  Ausuh-  However,  since  in  transition  metals 
the  on-site  energy  of  the  impurity  B  are  often  affected  the  host  A,  it  is  more 
accurate  to  find  the  impurity  potential  by  fitting  the  tight-binding  model  to  the 
density  functional  theory  for  a  small  unit  cell  such  as  A3B1.  Again,  the  error 
between  the  tight-binding  and  density-functional  eigenvalues  given  in  Equation 
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2.58  is  minimized.  Since  in  the  transition  metals  the  most  important  difference 
between  different  elements  is  the  d-orbital  on-site  energy,  all  the  parameters  except 
for  the  d-orbital  on-site  energy  of  the  impurity  can  be  set  as  the  parameter  of  the 
host  obtained  from  the  bulk  fit.  As  a  result,  the  fitting  for  A3B1  is  just  an  one- 
parameter  fit  for  each  spin,  where  the  only  varying  parameter  is  the  d-orbital 
on-site  energy. 

The  validity  of  the  parameters  are  then  confirmed  in  small  unit  cells  by  com- 
paring the  density  of  states  computed  by  the  tight-binding  theory  and  the  density- 
functional  theory.  For  example,  parameters  for  Fe  as  impurities  in  fee  Ni  host  are 
obtained  from  fits  to  fee  NisFe  band  structures.  As  a  check,  parameters  obtained 
from  the  fit  are  used  to  calculate  the  tight-binding  band  structure  of  fee  NisiFe, 
which  agrees  with  the  band  structure  obtained  from  the  density-functional  theory. 
This  indicates  that  the  model  works  well  at  least  when  impurity  concentration  is 
less  than  25%.  When  the  impurity  concentration  is  large,  or  when  the  difference 
in  hoping  between  the  host  and  the  impurity  is  large,  the  error  is  expected  to 
increase. 

The  other  parameter  in  interface  diffusion  is  the  density  of  impurity  atoms  in 
the  few  atomic  layers  next  to  the  impurity.  Unfortunately  this  parameter  depends 
the  material  and  methods  used  to  grow  the  sample  in  experiments.  Only  in  a  few 
transport  experiments  the  sample  is  characterized  well  enough  to  give  the  details 
of  the  interfaces.  Even  in  those  case  when  the  interfaces  are  well  characterized, 
such  as  those  by  Cyrille  et  al.  [6-7],  only  the  long  length  scale  («10  nm)  interface 
roughness  is  measured.  Within  our  model,  the  interface  roughness  at  this  length 
scale  should  be  considered  in  the  geometry  because  the  Fermi  wavelength  is  an 
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order  of  magnitude  smaller.  On  the  other  hand,  the  atomic  scale  interfaces  rough- 
ness is  related  to  the  density  of  impurity  atoms.  In  Chapter  4,  we  will  explain  how 
to  set  the  density  of  the  impurity  atoms. 

2.3.3  Impurity  Potential  for  Structural  Disorder 

The  impurity  potential  for  the  bulk  structural  disorder,  Awst,  can  be  obtained 
by  comparing  the  calculation  of  the  resistivity  and  the  resistivity  for  bulk  samples. 
Usually,  the  samples  in  heterostructures  are  not  as  clean  as  a  bulk  sample.  There- 
fore, we  use  the  room  temperature  resistivity  in  bulk  samples  as  a  guide  to  setting 
the  bulk  disorder  parameter. 

2.3.4  Extension  to  Magnetic  Alloys 

The  method  described  in  this  section  can  also  be  extended  to  model  transition 
magnet  alloys,  where  the  magnetic  moment  and  hence  the  band  structiure  vary  with 
the  composition.  In  these  materials,  the  localized  d  electrons  are  responsible  for 
the  magnetism.  We  can  assume  that  the  most  important  changes  due  to  alloying 
are  contained  in  the  on-site  energies,  i/(r,  r,  d,  d,  t,  t)  and  /f(r,  r,  d,  d,  |,  |),  of  the 
spin-up  (majority)  and  spin-down  (minority)  d  electrons.  In  other  words,  the  alloy 
is  assumed  to  have  site-independent  hoping  parameters  and  s  and  p  on-site  energies 
the  same  as  the  host  metal,  and  site-dependent  d  on-site  energies.  This  assumption 
will  be  justified  by  the  agreement  with  local  density  calculations  of  supercells  in 
different  impurity  concentrations. 

The  major  contributions  to  the  on-site  energies  come  from  the  atomic  core 
potential  plus  the  Coulomb  energy  due  to  the  opposite  spin.   In  the  itinerant 
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electron  model,  it  is  more  convenient  to  work  with  the  total  number  of  spin-up 
and  spin-down  d  electrons,  N}^  and  A'^^,  at  each  site  r.  Therefore,  we  rewrite  the 
parameters  in  the  form  of 

F(r,  r,  rf,    r,  T)   =   U^  +  U^Ni^d 

H{r,r,d,d,i,i)   =   U^  +  U^Nl„  (2.59) 

where  [7°  and  Up  are  the  on-site  parameter  and  effective  Coulomb  energy  per  pair 
at  site  r.  Both  A^J^  and  A/^f^  in  Equation  2.59  are  obtained  from  the  band  structure, 
which  in  turn  depends  on  .^^(r,  r,  d,  d,  j,  t)  and  H{r,  r,  d,  d,  |,  j).  Thus,  Equation 
2.59  has  to  be  solved  self-consistently  with  the  alloy  band  structure  calculated 
using  Hamiltonian  of  Equation  2.1.  The  alloy  band  structmre  is  computed  with 
the  coherent  potential  approximation  (CPA)  [49-50]  is  used,  which  is  summarized 
in  Appendix  B. 

All  parameters  for  the  host  atoms  can  be  found  from  fits  to  the  local  density 
calculations  of  the  pure  metal  [48].  The  only  two  parameters  left,  C/j^p  and  11^^^  of 
impm-ity  atoms,  are  obtained  from  fits  to  the  local  density  calculations  of  supercells 
composed  the  two  types  of  atoms. 

In  Figm-e  2.3,  we  plotted  the  Slater-Pauling  curve,  which  shows  the  magnetic 
moment  against  the  average  number  of  electrons,  for  3d  magnetic  alloys.  Both 
experimental  data  and  tight-binding  CPA  calculations  are  shown.  As  we  can  see, 
the  theory  agrees  not  only  with  the  main  curve,  but  also  the  NiCr  branch  as  well. 
The  NiCr  curve  is  very  sensitive  to  the  details  of  the  band  structure.  In  this 
simple  model,  only  two  parameters  in  addition  to  the  parameter  of  the  host  are 
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included;  therefore,  there  is  some  deviation  in  the  NiCr  curve.  The  agreement 
indicates  that  the  tight-binding  model  described  here  is  suitable  for  obtaining  the 
magnetic  properties  and  band  structure  properties  of  transition  magnetic  metal 
heterostructures. 
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Figm-e  2.3:  The  magnetic  moment  of  alloys  per  atom  versus  the  number  of  valence 
electrons  per  atom.  The  lines  are  the  experimental  data  of  FeV,  NiFe,  NiCu,  and 
NiCr  alloys  as  labeled.  The  symbols  are  from  the  CPA  calculations  of  FeCr  (tri- 
angles), NiCr  (circles),  NiCu  (squares),  and  NiCr  (crosses).  The  CPA  calculation 
agrees  well  with  the  experiment  in  the  main  branch  (NiFe  and  NiCu),  and  agrees 
reasonably  well  with  the  NiCr  branch.  The  dip  in  the  experiment  curve  in  NiFe 
corresponds  to  the  bcc-fcc  phase  transition. 


CHAPTER  3 
COMPUTATION  PROCEDURE 


In  the  last  chapter,  we  discussed  the  formal  aspects  of  modeUng  magnetic  nanos- 
tructnres  and  calculating  the  conductivity.  We  also  discussed  how  to  obtain  pa- 
rameters so  that  the  models  give  reaUstic  results  that  can  be  used  to  compare  with 
experiments. 

However,  if  one  writes  computer  codes  directly  from  these  equations,  the  re- 
sulting codes  will  not  be  fast  enough  for  large  scale  studies.  In  particular,  many 
of  the  equations  involve  0{N^)  operations  such  as  the  matrix  multiplication  of 
NM  X  NM  dense  matrix,  where  A*"  is  the  niunber  of  atoms  in  a  unit  cell,  and  M 
is  the  number  of  orbitals  per  site.  In  order  to  speed  up  the  calculation  to  tackle 
the  tasks  such  as  optimization  and  studying  large  systems,  we  have  to  reformulate 
the  problem  with  special  consideration  about  the  computation  speed. 

In  this  chapter,  we  discuss  the  computation  procedure  of  the  conductivity  cal- 
culation. In  the  first  section,  we  discuss  the  current-perpendicular-to-plane  (CPP) 
conductivity  computation,  which  demonstrates  how  to  find  the  impurity  averaged 
Green's  function,  solve  for  the  vertex  correction,  and  compute  the  conductivity. 
Although  we  use  the  CPP  geometry  as  an  example,  the  application  to  general 
geometry  are  similar.  Most  importantly,  the  tools  that  are  developed  here  can  be 
used  in  more  general  cases  and  the  size-dependence  of  the  computational  cost  are 
the  same.  As  we  will  show,  the  size-dependence  for  computing  the  self-energy,  the 
vertex  correction,  and  the  conductivity  are  0{N),  OiN"^),  and  0{N),  respectively. 
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This  is  an  0{N)  saving  compared  to  total  computation  costs  in  the  last  Chapter. 
As  a  result,  we  can  study  more  systems  with  a  larger  size.  As  we  will  see  in  Chapter 
5,  the  speed  of  the  code  is  crucial  to  optimization. 

In  the  second  section,  we  will  discuss  the  current-in-plane  (CIP)  geometry, 
which  is  also  commonly  studied  in  experiments.  We  will  show  that  since  the 
Hamiltonian  has  translational  symmetry  in  the  in-plane  direction,  the  vertex  cor- 
rection is  identical  to  zero.  Therefore  the  CIP  conductivity  is  easier  to  computer. 
However,  the  size-dependence  of  the  overall  computation  remains  0{N^). 

In  the  third  section,  we  briefly  discuss  the  calculation  of  the  conductivity  tensor 
and  the  current-at-an-angle-to-plane  (CAP)  geometry,  which  will  only  require  tools 
developed  in  the  first  two  sections. 

In  the  last  section  of  the  chapter,  we  explain  how  we  implement  the  algorithms 
discussed  in  this  chapter.  We  shall  use  some  new  C-|— I-  matrix  hbraries,  which 
will  call  the  FORTRAN  libraries  BLAS  (Basic  Linear  Algebra  Subprograms)  [51] 
and  LAPACK  [52]  for  the  actual  computation.  We  will  also  use  an  in-house  block- 
matrix  library  to  implement  most  algorithms,  and  explain  the  advantage  of  using 
this  new  tools  for  coding.  - 

3.1    Current-Perpendicular-to-Plane  Conductivity 

We  discuss  the  computational  aspects  of  the  current-perpendicular-to-plane 
(CPP)  conductivity  calculation  in  this  section.  Let  us  consider  a  periodic  multi- 
layer system  that  has  N  mono-atomic  layers  in  each  period  and  can  be  generated 
by  a  unit  cell  of  A''  atoms.  We  label  the  growth  direction  as  z  and  the  in-plane 
directions  as  x  and  y.  Since  in  the  imit  cell  each  atom  belongs  to  a  different  mono- 
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layer,  we  label  the  an  atom  in  the  unit  cell  as  z  =  1,2, . . . ,  N.  An  example  of  this  is 
a  multilayer  grown  in  the  fcc(lll)  direction.  Each  elements  of  a  matrix  is  an  m  x  m 
block  matrix,  where  m  is  the  nimiber  of  orbitals  per  site.  The  current  is  in  the 
^-direction.  The  equations  in  this  section  are  intended  for  direct  implementation, 
therefore,  the  computational  cost  of  major  equations  are  given. 

3.1.1    Impurity  Averaged  Green's  Function 

Using  the  model  discussed  in  Chapter  2,  the  self-energy  is  given  by  : 

S«(ri,r2)   -  t;(ri)C;^(ri,ri)(J(ri,r2) 

(3.1) 

We  only  need  to  calculated  the  real  space  function  for  sites  within  the  unit  cell. 
Since  S  is  diagonal,  it  can  be  stored  in  a  vector  and  the  computation  cost  is  N. 
In  real  space,  the  Green's  function  is  given  by 

G^={uj-H-  E^)-^  (3.2) 

However,  since  a  periodic  system  is  also  infinite,  we  have  to  calculate  the  Green's 
function  for  the  imit  cell  in  the  Fourier  space: 


(3.3) 
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where  k  is  a  wave- vector  in  the  first  Brillouin  zone,  i/k  and      are  the  Hamiltonian 
and  self-energy  in  the  k-space.  Since  the  self-energy  is  on-site,  we  have  = 
in  the  unit  cell. 

The  Hamiltonian  in  k-space  can  be  found  by  i/k  =  Hr -^^(0,  R)exp(— ik  •  R), 
where  R  is  the  position  vector  of  neighbor  unit  cells.  In  the  following,  we  consider 
up  to  second  neighbor  hoping  in  real  space.  In  the  fcc(lll)  lattice,  up  to  the 
second  nearest  neighbors  are  included  in  the  three  monolayers  described  by  ^  = 
0,  ±1.  Therefore,  when  we  work  in  the  k-space  representation,  the  Hamiltonian  is 
tridiagonal  with  periodic  boundary  condition.  In  Appendix  A,  we  explain  how  the 
inverse  can  be  found  efficiently.  Since  ffk  is  tridiagonal  with  periodic  boundary 
condition  and  E  is  diagonal,  the  cost  for  inverse  is  0{N'^).  However,  since  only 
diagonal  elements  of  G  are  needed  to  find  the  self-energy,  the  cost  can  be  cut  to 
0{N)  in  this  step. 

The  real  space  Green's  function  within  the  unit  cell  is  then  given  by 

G^  =  ^EGk  (3.4) 

■(Vk  k 

The  advanced  self-energy  and  Green's  function  can  be  foimd  by  taking  the  hermi- 
tian  conjugates  of  the  results.  The  Equations  3.1  and  3.4  are  solved  self-consistently 
by  iteration. 

If  the  system  is  finite  instead,  then  we  can  work  in  the  real  space  and  the  inverse 
algorithm  for  fixed  boundary  condition  in  Appendix  A  can  be  used;  however,  one 
has  to  consider  the  lead  with  special  care. 
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3.1.2    Vertex  Correction 

To  find  the  conductivity,  we  have  to  compute  not  only  the  diagonal  of  the 
Green's  function,  but  also  the  full  Green's  function.  In  Appendix  A,  we  explain  how 
to  compute  the  Green's  function  in  0{N'^)  for  both  periodic  and  fixed  boundary 
conditions. 

In  the  k-space,  the  current  density  operator  at  site  z  in  the  unit  cell  is  given 

byi 

Jzk{zi,  Z2)   =   {-ia3/V,)Hi,(zi,  Z2)      a  zi=  z,Z2  =  z  +  1 
^   {ia3/Vz)Hk{zi,  Z2)      a  zi  =  z  +  l,Z2-^  z 
=   0  otherwise, 

(3.5) 

where  03  is  the  monolayer  distance.  The  total  current  is 

Jk  =  E'^^kK,  (3.6) 

z 

which  has  2A''  elements. 

As  discussed  in  Subsection  2.2.3,  both  sets  of  equations.  Equations  2.32-2.37 
and  2.38-2.43,  can  be  used  to  compute  the  conductivity.  When  we  compute  the 
vertex  correction,  both  Equations  2.32  and  2.38  have  the  same  computational  costs. 
To  see  this,  we  notice  from  the  definition  of  K^^  that  it  is  non-local.  Therefore, 
both  K^^  and  K^"^  have  the  same  number  of  elements  and  are  computed  in  the 

^Here  the  coordinates  are  cyclic,  for  example,  z  =  N  +1  refers  to  the  same  site  as  2  =  1. 
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same  number  of  steps.  Since  is  local,  the  cost  for  computing  lix[GTGJz]  in 
Equation  2.43  is  0{N),  while  the  cost  of  Tr[GJGr,]  in  Equation  2.37  is  0{N'^). 

The  vertex  correction  is  given  by  the  Fourier  transform  of  the  Dyson's 
equation  (Equation  2.36): 


k  24,25 

+  E^(^i)^k  (^1,  z,)K^\z,,  z,)G^{z,,  z,) 

Z4 

=^    ^i^l)  E  E[<^k  (^1,  2:4)4(24,  Z4  +  1)G^{Z4  +  1,  21) 
k  24 

+G^izi,ZA  +  1)4(24  +  1,  z)G^iz4, 21)] 
+^(^i)E[G'k(^i' ^4)^^^(24,24)^1^(24, 21)]  (3.7) 
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Again,  K^^  =  K^"^  in  the  unit  cell  because  K^^  is  diagonal.  The  first  and  second 
term  cost  0{N'^),  however,  they  are  independent  of  the  iteration.  In  the  last  term, 
for  each  21,  we  need  to  do  a  sum  over  24,  thus,  each  iteration  in  this  loop  costs 
0{N^).  Since  this  step  is  inside  the  iteration  loop,  it  is  the  most  time-consuming 
part  of  the  whole  calculation.  The  cost  for  computing  G^  to  start  the  loop  is  also 
0{N^). 

3.1.3  Conductivity 

The  conductivity  is  given  by  Fourier  transforming  Equation  2.43  as 

=   JzET^  [24kGi^(4  +  K^^)G^  -  4kG^4Gfk  -  J.uGUCt]  -(3.8) 
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As  explained  in  Chapter  2,  is  independent  of  2,  and  hence  we  can  choose  any  z 
the  evaluate  the  conductivity. 

One  should  note  that  it  only  take  0{'N)  to  evaluate  the  traces.  For  example, 


the  first  two  terms  can  be  written  as 


Tr 


+{G^J,y,G^){z',z'  +  l)Mz'  +  l,z'),  '   '  (3.9) 


and 


V.  J, 


Tr 


GUi^G^K^^]   =  Y.iGtJz^^G^){z\z')K^\z\z'l  (3.10) 


where 


{GtJ.^.G^'\z'-^\,z')  = 


{GtJ,^GZ'\z\z'  +  l)  = 


{GUyG^^^)iz',z')  = 


G^{z'  +  1,  z)J,y,{z,  z  +  l)G^^^{z  +  1,  z') 
+G^{z'  +  1,2  +  1)  J,k(^  +  1,  z)G^^^(z,  z'), 
G^{z',  z)J,^{z,  z  +  1)GZ'^{z  +  1,2'  +  1) 

+Gi^(2',  2  +  1)  J,k(2  +  1,  2)Gf/^(2,  2'  +  1), 

Gt{z\z)J,y,{z,z  +  l)G^/''{z  +  l,z') 
+G^{z',z  +  l)J,^{z  +  \,z)G^"'{z,z').  ' 


(3.11) 


Note  that  there  is  only  one  sum  in  each  trace,  therefore,  the  cost  for  each  trace  is 
0{N).  Further  simplifications  are  possible  by  using  the  symmetry  of  the  matrices. 


for  example, 
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Tr[(G^J,G^)(^',  z'  +  l)J{z'  +  =  Tr{[J{z',  z'  +         J,G«)(z'  +  1,  z')]^ 

=   (Tr{(G^ J,G«)(/  +  1,  z')J{z\  z'  +  1)})*. 

(3.12) 

Therefore, 

J,G^  J]  =  2  ^  a^KG^  J,G'^)(z'  +  1,  z')       z'  +  1)}.  (3.13) 

3.1.4    Diagonal  Impurity  Potential 

Most  of  the  discussions  in  this  work  allow  the  impurity  potential  to  be  off- 
diagonal  in  the  orbital  indices.  However,  for  the  diagonal  impurity  potential, 
v{z){m,n)  =  v{z){m)5{m,n),  further  simplification  that  reduces  the  cost  of  the 
second  loop  to  negligible  is  are  possible.  In  this  case,  since  K^"^  is  also  diago- 
nal in  band  index,  we  denote  K^^{z){m)  =  K^^{z,  z){m,m).  In  the  following, 
we  will  show  that  the  last  term  in  Equation  3.7  can  be  written  in  the  form  of  a 
matrix- vector  multiplication.  Hence  we  can  group  the  K^^  terms  on  both  sides  of 
Equation  3.7  and  write  this  Dyson's  equation  in  the  form  of 

il-q)K^^=p,  (3.14) 

where  g  is  a  matrix  and  p  is  a  vector.  Equation  3.14  is  a  set  of  hnear  equations 
and  can  be  solved  without  iterations.  To  derive  this  equation,  we  notice  that  for 
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each  zi, 

K'^{z,){m)   =  p{z,){m) 


and 


-^^^^  Y:  E  G§{z„Z2){m,  n)K{z2){n)Gt{z„  z,){n,  m) 


k  22 


=  p{z,){m)  +  Y.q{z,,Z2){m,n)K''\z2){nl  (3.15) 


22 


where 


^^^^  E  E[G^(^i>  ^2)  Jk(^2,  ^2  +        (^2  +  1,  ^l) 


k  22 


+Gk(-2l,  22  +  1)4(22  +  1,  22)Gk  (22,  2l)](m,  m) 

^^^^^EE[Gk(^i>^2)Jk(22,22  +  l)Gi^(22  +  l,2i)](m,m)  +  h.C. 

k  22 

E  E  »{[Gl^(2l,  22)4(22,  22  +  1)G^^(^2  +  1, 2i)](m,  m)}. 

(3.16) 


g(zi,^2)(m,n)   =  ^  g«(zi,  Z2)(m,  n)Gi^(^2, 2i)(n,  m) 

^  k 

y(2i)(m) 


El<^k(^i'^2)(m,n)|l  (3.17) 
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Since  the  second  term  of  p  is  the  conjugate  of  the  first  term,  the  only  term  in  p  we 
need  to  compute  is  the  first  term.  It  can  be  calculated  as  follows  : 


[G^{zu  Z2)Mz2,  Z2  +  1)G^{Z2  +  1,  Zi)]{m,  m) 
=  E[<^k  (2i>  Z2)Mz2,  ^2  +         m')G^iz2  +  1,  zi){m',  m), 

m' 

(3.18) 

which  does  not  involve  any  sum  over  positions.  Since  q  is  symmetric,  we  only  need 
to  calculate  half  of  the  element,  and  we  can  use  the  symmetric  matrix  solver  for 
this  problem.  These  two  cost  saving  steps  are  important  because  they  are  0{N^) 
calculations. 

When  impiurity  potential  is  not  diagonal  in  the  orbital  indices,  K^"^  and  G  do 
not  commute  and  it  is  very  costly  to  solve  the  linear  equations,  therefore  K^"^ 
is  solved  by  iteration.  However,  for  the  impurity  potential  used  here,  it  is  much 
faster  to  Equation  3.14  directly,  which  is  a  set  of  real  and  symmetric  linear  equa- 
tions and  can  be  solved  with  a  cost  of  The  cost  for  computing  p  and  q 
are  proportional  to  NkN'^M^  and  NkN'^M'^,  where  Nk  is  the  number  of  k-points. 
Therefore,  most  time  is  spent  in  evaluating  p. 

3.1.5    Total  Computational  Cost 

The  computational  cost  is  important  to  the  large  scale  we  discuss  in  the  follow- 
ing Chapters,  therefore  it  is  given  here.  If  A^^,  are  the  number  of  iterations  in  the 
self-energy  loop  respectively,  and  Nk  is  the  number  of  k-points,  then  the  total  cost 
for  a  conductivity  calculation  is  roughly  {4M^)Nk{6N^  +  {l9Ns  +  27) N)  floating 


point  operations  (flop),  where  the  factor  4M^  is  the  cost  for  a  complex  9x9  block 
matrix  multipUcation.  As  an  estimate,  we  take  Nk  ^  1700,  M  =  9,  A''  =  50, 
Ng  —  10,  the  total  cost  is  on  the  order  of  10^^  flop.  On  a  machine  capable  of 
about  10^  flop  per  second  (1  GFlops),  this  should  take  a  few  minutes.  In  practice, 
we  performed  the  calculations  a  computer  cluster  with  ten  866MHz  Pentium  III 
CPUs,  it  takes  about  5  to  10  minutes  for  a  calculation  when  jV  =  50. 

The  algorithms  presented  in  this  chapter  are  CPU  intensive.  The  memory 
requirement  is  0{N^).  For  the  computers  used  in  this  work,  unless  N  is  larger 
than  about  500,  one  does  not  need  to  rearrange  terms  to  optimize  the  memory 
usage. 

3.2    Current-In-Plane  Conductivity 

The  current-in-plane  (CIP)  conductivity  is  a  special  case  where  the  Hamiltonian 
has  a  translational  symmetry  in  the  direction  of  the  current.  Let  us  consider  the 
same  fcc(lll)  lattice  discussed  in  the  last  section  and  denote  the  direction  of  the 
current  as  x.  The  problem  is  greatly  simplified  because  the  vertex  correction  AT^^ 
is  zero.  This  result  can  be  shown  by  studying  the  first  correction  term  of  Equation 
2.40: 

(^G^7G^)(x,x)   =  v{x)[J2G''{x,x  +  j)Jix+j,x  +  j  +  l)G^{x  +  j  +  l,x) 

j 

+  ^         x  +  j  +  l)Jix  +  j  +  l,x+  j)G\x  +  3,  x). 
j 

(3.19) 
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In  the  real  space  representation,  the  Hamiltonian  is  real  and  has  translational 
symmetry,  hence  and  are  symmetric,  and  J  is  antisymmetric.  The  second 
sum  above  can  be  rewritten  as 

-  J2  G^i^  +  i  +  1,  r)J{x  +  j,x  +  j  +  l)G^{x,  X  +  j) 
J 

=       G^^i^^  ^  -  i  -  ^Ui^  -j-hx-  3)G\x  -  i,  x),  (3.20) 

i 

where  in  the  last  line  we  have  iised  the  translational  symmetry  separately  for  G^, 
J,  and  G^.  By  substitute  /  =  j  +  1,  we  can  see  the  two  sums  in  Equation  3.19 
cancel  exactly.  As  a  result,  all  the  correction  terms  to  F^^  in  Equation  2.40  are 
also  zero. 

Since  the  vertex  correction  is  zero,  we  can  compute  the  conductivity  corre- 
sponding to  the  open  bubble  by 

c^a;  =  —  ^Tr[2JkG'j^Jj;kG^  -  JkG^JxkG^  -  JkG^JxkG^].  (3.21) 

Since  J^k  is  a  tridiagonal  matrix  in  the  unit  cell,  it  has  SN  elements.  Thus,  the 
computational  cost  for  each  trace  is  0(A''^). 

3.3    Conductivity  Tensor 

As  explained  in  Chapter  2,  the  current  density  operator  is  not  obvious  to  con- 
struct. It  is  especially  difficult  when  the  direction  of  cmrrent  is  at  an  angle  to  the 
crystal  axes.  Therefore,  in  computing  the  cmrent  for  systems  such  as  the  current- 


at-an-angle-toplane  (CAP)  system,  it  is  easier  to  compute  the  conductivity  tensor 
and  use  the  tensor  equation  J  =  crE. 

In  a  cubic  bulk  system,  a°'^  =  Sa^cr^^-  However,  in  nanostructures,  the  sym- 
metry is  usually  broken,  therefore  a°'^  are  in  general  different.  The  techniques  for 
computing  current  with  vertex  correction  is  useful  for  computing  the  conductivity 
tensor. 

3.4  Coding 

As  shown  in  the  above,  it  is  natural  to  formulate  our  problem  in  block  matrices 
and  hence  nattual  to  write  our  code  using  block  matrix  object  classes.  We  have 
to  implement  the  block  matrix  algorithms  such  as  inverse  and  products  discussed 
in  previous  sections  to  take  advantage  of  the  sparsity  of  the  matrices.  Therefore, 
it  is  not  enough  to  just  have  a  block  matrix  class  and  common  algorithms.  The 
block  matrix  class  should  help  us  develop  new  block  matrix  algorithm  efficiently. 

Our  goal  is  to  translate  mathematical  equations  involving  block  matrices  into 
programming  lines  directly,  while  having  performance  close  to  Fortran  codes. 

In  this  section,  we  discuss  how  this  is  achieved  with  C+-I-.  First  we  use  the 
Template  Numerical  Toolkit  package  to  handle  the  operations  of  all  dense  matrices. 
To  increase  the  speed,  we  link  critical  subroutines,  for  example,  the  matrix  mul- 
tiplication or  inverse  operation  for  dense  blocks,  to  the  Fortran/ Assembler  coded 
native  binaries  BLAS  [51]  and  LAPACK  [52].  Next  we  develop  a  block  matrix  class 
that  can  use  the  TNT  matrix  class  as  elements.  Using  the  new  block  matrix  class, 
we  develop  codes  to  implement  the  algorithms  described  by  previous  sections. 
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In  addition  to  these  steps,  our  code  is  parallelized  with  PVM  (Parallel  Virtual 
Machine)  at  the  level  of  k-points  to  use  the  all  available  computing  resources  of  a 
heterogeneous  network.  The  resulting  code  should  scale  well  with  the  number  of 
processors. 

3.4.1    Template  Numerical  Toolkit  (TNT) 

The  template  [53]  feature  in  C++  language  allows  writing  a  single  code  for 
.  :         object  classes  or  algorithms  for  different  storage  types.  For  example,  the  subroutine 

template  <class  T> 
T  f  (T  x) 

iVri ,  return  x*x: 

-  } 

■  defines  a  function  that  can  be  called  by  f  (x)  whether  the  data  type  of  x  is  int 

-  •  or  complex<  double  >.  The  compiler  generates  suitable  binaries  to  match  for  the 

'  datatype. 

When  we  started  developing  the  first  version  our  code,  the  ANSI  C++  standard 
was  not  formally  established  yet.  However,  there  were  two  attractive  experimental 
C++  packages  that  handle  matrices  well  by  making  use  of  the  template  featiues: 
the  Template  Numerical  Toolkit  (TNT)  by  R.  Pozo  [54-55]  and  Blitz++  by  T. 
Veldhuizen  [56-57].  Although  Blitz++  has  great  potential,  it  used  advanced  fea- 
tures such  as  expression  templates,  which  were  not  supported  by  the  compiler  on 
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our  platforms  at  the  time  we  started  coding.^  On  the  other  hand,  the  TNT  pack- 
age could  be  compiled  and  linked  with  Fortran  programs  on  our  major  platforms: 
Solaris  and  Linux.  Consequently,  we  chose  the  Template  Numerical  Toolkit. 

The  TNT  object  we  use  most  are  the  Fortran_Matrix  object,  which  stores 
a  matrix  using  Fortran  (column)  convention,  and  the  Vector  object.  They  allow 
Fortran  (one-based)  indices  such  as  A (2,3) .  For  example,  the  following  lines  define 
a  2  X  2  matrix  A,  vectors  x  and  b. 

Fortran_Matrix<  int  >  A(2,2,  "12 

14"); 

Vector<  int  >  x(2,  "  12"),  b(2) ; 
b(l)=3;  b(2)=4; 

The  package  also  takes  advantage  of  the  function/operator  overload  feature  of 
C-I-+  to  make  the  matrix  and  vector  operations  intuitive.  Basic  operations  like 
=,  +,  *,  and  «  are  overloaded.  Continued  from  the  above  example,  the  next  two 
lines  computes  and  prints  out  the  vector  c  =  SA^x  +  b. 

Vector<  int  >  c=3*A*A*x+b; 
cout  «  c  «  endl; 

Notice  that  the  compiler  links  the  three  multiplications  above  to  scalar-matrix  mul- 
tiplication, matrix-matrix  multiplication,  and  matrix- vector  multiplication  subrou- 
tines in  order.  Without  operator  overloading,  the  line  that  computes  c  would  look 
like 

Vector<  int  >  c(vecadd(vecmult(matmult(scalmult(3,A) ,A) ,x) ,b)) ; 

2 According  to  the  Blitz+-|-  website  [57],  Blitz++  can  now  be  compiled  by  compilers  like  gcc 
and  KAI  C++. 
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3.4.2  Fortran  performance  with  BLAS  and  LAPACK 

To  achieve  the  performance  close  to  that  of  a  Fortran  code,  we  overload  the 
most  often  used  and  most  numerically  intense  subroutines  by  linking  them  to 
BLAS  and  LAPACK.  For  example,  we  overload  the  matrix  multiplications  A*B  by 
the  BLAS  matrix  multiplication  subroutines  xGEMM,  where  x  designates  the  data 
type  of  the  matrices.  This  is  done  efficiently  because  the  Fortraa_Matrix  object 
stores  matrix  element  using  Fortran  convention. 

To  make  our  code  easy  to  use  and  understand,  we  name  commonly  used  func- 
tions similar  to  Matlab  commands.  For  example,  to  find  the  inverse  of  a  ma- 
trix, we  use  B=inv(A),  which  is  linked  to  the  LAPACK  subroutines  xGETRF  (LU- 
factorization)  followed  by  xGETRI  (inverse  of  an  LU-factored  matrix). 

Since  owe  code  spends  most  of  the  time  doing  matrix  operations  implemented 
with  Fortran/ Assembler  native  binaries,  the  overall  performance  of  our  code  is 
close  to  a  pure  Fortran  code.  On  the  other  hand,  the  codes  are  easy  to  write  and 
understand. 

3.4.3  BlockMatrix  class 

We  design  a  BlockMatrix  class  to  implement  mathematical  algorithms  as  di- 
rectly as  possible.  Without  such  an  object  class,  block  matrix  algorithms  are 
difficult  to  code  and  debug.  In  the  following,  we  explain  which  featines  we  want  in 
a  BlockMatrix  class,  and  how  these  features  enable  us  to  implement  comphcated 
block  matrix  equations  without  much  effort. 

We  want  to  design  a  BlockMatrix  class  as  a  matrix  Fortran_Matrix  objects. 
Using  the  Fortran_Matrix  object,  the  elements  (blocks)  of  our  BlockMatrix  class 
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will  be  intuitive  and  fast  to  operate.  The  storage  should  be  efficient,  which  means 
that  block  matrices  of  different  sparsity  have  different  storage  arrangements  and 
different  index  functions.  To  do  this,  we  define  several  subclasses  such  as  the 
Full_BlockMatrix,  PeriodicBand-BlockMatrix,  and  the  Sparse_BlockMatrix. 
The  index  functions  are  defined  within  the  subclasses  to  overwrite  the  virtual 
index  functions  defined  in  the  parent  class. 

The  object  should  be  intuitive  to  use.  For  example,  one  may  obtain  the  current 
operator  from  a  tight-binding  Hamiltonian  with  the  following  code: 

for  (int  r=l;  r<=N-l;  r++) 
{ 

J(r+l,r)=i*H(r+l,r); 
J(r,r+l)=dagger(J(r+l,r)) ; 

} 

Here,  H(j  ,  j+1)  is  a  Fortran_Matrix  itself,  and  also  an  element  of  the  BlockMatrix 
H.  In  addition,  we  can  access  elements  of  the  matrix  H(j ,  by  an  expression 
like  H(j,j+l)(ml,ni2). 

These  featm-es  allow  us  to  implement  block  matrix  algorithms  efficiently  and 
more  importantly,  debug  in  a  shorter  time. 

In  the  following,  we  show  how  the  of  Godfrin  inverse  method  is  implemented 
using  the  BlockMatrix  class.  We  want  to  find  G  =  where  M  is  a  tridiagonal 
block  matrix.  These  matrices  are  stored  as  PeriodicBand-BlockMatrix,  which 
is  a  subclass  of  BlockMatrix  that  stores  a  banded  block  matrix  with  periodic 
boundary  conditions.  For  example,  M  is  declared  by 
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PeriodicBand_BlockMatrix<  Fortran_Matrix<  complex<  double  >  >  > 
M(N,  KL,  KD,  KU,  p,  q) ; 

where  N  is  the  number  of  blocks,  KL  =  KD  =  KU  =  1  are  the  number  of  bands  in  the 
lower  triangle,  diagonal,  and  upper  triangle  or  the  matrix,  and  the  block  size  is 
p  X  q. 

The  Godfrin  method  shown  in  Equations  A. 21  and  A.  11  finds  the  inverse  of  a 
tridiagonal  block  matrix  M  in  0{N'^).  To  implement  this  algorithm,  we  first  declare 
the  subroutine,  with  input  matrix  M,  and  output  matrix  G  and  other  variables  used 
as  temporary  memory.  Values  such  as  the  matrix  size  and  block  size  are  initialized, 
as  shown  below.^ 

template  <class  T> 

GodfrinCconst  PeriodicBand_BlockMatrix<  T  >  &M, 
PeriodicBaiid_BlockMatrix<  T  >  &G, 
Full-BlockVector<  T  >  &X,  Full_BlockVector<  T  >  &Y, 
Full_BlockVector<  T  >  &C,  Full-BlockVector<  T  >  &D) 

{ 

int  N=M.num_cols() ; 
int  p=M.block_h() ; 
int  q=M.block_w() ; 
T  W(p,q); 

Here,  class  T  is  the  class  of  an  element  of  a  block  matrix,  which  is  Fortran-Matrix< 
complex<  double  >  >   in  this  example. 

^For  clarity,  the  bound  checking  codes  are  not  shown.  Also,  the  actual  code  is  slightly  modified 
to  output  either  the  full  matrix  or  selected  bands  of  G  depending  on  the  input. 


The  next  code  segment  is  directly  translated  from  Equation  A. 21. 


67 


X(N)=0.0; 

for  (int  j=N;  j>=2;  j~) 
{ 

W=inv(X(j)-M(j,j)); 

C(j)=W*M(j,j-l); 

X(j-l)=neg(M(j-lJ)*C(j)); 


Y(1)=0.0; 

for  (int  j=l;  j<=N-l;  j++) 
{ 

W=inv(Y(j)-M(j,j)); 

D(j)=W*M(j.j+l); 

Y(j+l)=neg(M(j+l.j)*D(j)); 

} 

Since  in  Equation  A. 21,  Wj  ans  Zj  are  not  used  later,  the  memory  of  W  is  reused. 

Notice  that  variables  like  M(j-l,j)  and  W  are  themselves  matrices.  Although 
M  is  tridiagonal,  there  is  no  need  to  worry  about  where  M(j-l,j)  is  stored  when 
writing  this  code.  The  details  of  the  storage  are  handled  in  the  index  function  of 
the  BlockMatrix  object  class  according  to  the  sparsity  of  the  M.  What  is  left  is 
straight  forward  translation  of  equations.  This  code  segment  should  be  transparent 
even  to  people  without  C++  experience,  and  it  is  therefore  easy  to  debug. 


The  last  part  of  the  subroutine  computes  the  Green's  function  using  Equation 
A.  11.  Again,  this  is  a  direct  translation. 

for  (int  j=l;  j<=N;  j++) 
{ 

G(j.j)=inv(M(j.j)-X(j)-Y(j)); 

for  (int  jl=j+l;  jl<=N;  G(jl, j)=C(jl)*G(jl-l, j) ; 

for  (int  jl=j-l;  jl>=l;  jl--)  G(jl, j)=D(jl)*G(jl+l, j) ; 

} 

} 

The  time  spend  on  designing  the  BlockMatrix  class  should  be  less  then  the  time 
we  saved  in  coding  and  debugging  other  part  of  the  code.  Since  the  BlockMatrix 
object  classes  developed  in  this  project  are  reusable,  future  implementations  of 
mathematical  algorithms  involving  block  matrices  should  be  as  easy  as  the  example 
shown  above. 


CHAPTER  4 

APPLICATION  TO  GIANT  MAGNETORESISTANCE 


In  the  previous  chapter,  we  explained  how  to  compute  the  transport  properties 
in  magnetic  nanostructures.  In  this  chapter,  we  apply  the  method  to  compute 
the  current-perpendicular-to-piane  giant  magnetoresistance  (CPP-GMR)  [4-14]  in 
magnetic  multilayers. 

As  explained  before,  a  good  theory  has  to  agree  with  both  the  conductivity  and 
the  GMR.  To  date,  we  have  found  two  sets  of  experiments  on  the  Fe/Cr  CPP-GMR 
with  both  published  resistivity  and  GMR  data,  one  is  by  Gijs  et  al.  [4]  and  the 
other  is  by  Cyrille  et  al.  [6-7].  In  the  following,  we  first  review  these  experiments. 
The  we  try  to  see  if  we  can  use  theoretical  models  to  reproduce  the  experimental 
data  for  reasonable  values  of  the  parameters. 

In  order  to  understand  the  trends  shown  in  the  experiments,  the  magnitude  of 
the  computed  resistivity  be  similar  to  the  values  found  in  the  experiments.  The 
changes  in  resistivity  and  GMR  as  a  function  of  different  parameters  are  much 
smaller  in  scale,  they  should  be  studied  after  the  magnitude  of  the  resistivity  is 
understood. 

Using  the  model  in  Chapter  2,  we  are  able  to  calculate  both  the  resistivity 
and  GMR  values  which  axe  close  to  Cyrille's  data  within  reasonable  parameters. 
However,  the  resistivity  data  of  Gijs  seem  to  be  too  small  to  fit  in  our  theory  unless 
very  small  disorder  is  used.  We  believe  the  required  disorder  is  too  small  to  be 
physical. 
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4.1    Experiment  of  Gijs  and  Collaborators 

The  Fe/Cr  superlattice  in  the  experiment  Gijs  et  al.  [4-5]  are  micro-structured 
pillars,  shown  schematically  in  Figure  1.2.  The  pillar  in  the  middle  in  Figure  1.2(d) 
is  the  sample  to  be  measm-ed.  The  samples  sizes  and  measured  conductivity  and 
OMR  are  summarized  in  Table  1. 


Table  4.1:  Summary  of  Gijs  [Fe(30  A)/Cr(tcr)]ioo  samples 


Cr  thickness  ^cr 

A 

Area 

GMR 

Psat(r  =  0) 

10 

90 

1.08 

11.0 

28 

20 

0.08 

24.3 

40 

6 

0.04 

25.2 

The  data  are  presented  in  terms  of  GMR  =  Ro/Rsm  -  1,  shown  in  Figure  4.1, 
and  the  quantity  \J{Ro-  Rsa.t)RoA'^,  shown  in  4.1.  Using  these  two  quantities  and 
the  areas  given  in  table  4.1,  we  can  obtain  the  resistance  Rq  and  i?sat,  and  hence 
po  and  psat.  At  low  temperatiue,  the  resistivity  of  the  samples  is  about  11  fiUcm 
for  Cr  layer  thickness  tcr  =  lOA,  and  is  increased  to  about  25  nilcm  for  ^cr  =  28A 
and  40  A.  The  GMR  of  the  sample  is  108%  for  ^cr  =  lOA  and  drops  to  4-8%  for 
tcr  =  28A  and  40  A. 


71 


As  we  shall  see,  the  resistivity  for  both  saturation  field  and  zero  field  are  very 
important  to  comparing  with  the  theory. 

120  I  , 


T(K} 

Figure  4.1:  Temperature  dependence  of  CPP  and  CIP  GMR  of  pillar  structures 
by  Gijs  et  al.  [4]. 


4.2    Experiment  of  Cyrille  and  Collaborators 


Cyrille  et  al.  [6-7]  studied  the  effect  of  the  CPP-GMR  as  the  number  of  bilayers 
N  of  the  Fe/Cr  superlattice  increases.  They  found  that  the  GMR  increases  as  A'' 
increases. 
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Figure  4.2:  Temperatiire  dependence  ^J{Ro-  /?sat)-Ro^^  for  in  the  data  of  Gijs  et 
al.  [4],  where  A  is  the  area  of  the  sample.  This  quantity  is  used  in  the  Valet- Pert 
theory  [16]. 

Figures  4.2  show  the  resistivity  and  the  GMR  of  a  series  of  [Fe(30  A)/Cr(12 
A)]jv  superlattices,  where  A'^  ranges  from  10  to  60.  In  this  figure,  the  resistivity 
for  both  saturation  field,  Psat,  and  zero  field,  po,  is  plotted  as  a  function  of  the 
number  of  bilayers  A'^.  The  magnitude  of  the  GPP  resistivity  is  about  35  fiQcm 
for  for  saturation  field,  and  is  about  40-50  fxi^cm  for  the  zero  field. 

The  experiments  also  show  that  psat  is  roughly  independent  of  A''  and  po  in- 
creases as  A''  increases.  This  trend  is  not  very  clear  on  Figiu-e  4.2(a),  probably  due 
to  sample  variations.  However,  since  po  and  psat  are  correlated,  the  trend  is  more 
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Figure  4.3:  CPP  Resistivities  and  GMR  in  FeCr  superlattice  as  a  function  of  num- 
ber of  bilayers  measured  by  Cyrille  et  al.  [6-7].  In  (a),  the  zero  field  resistivity  po 
increases  as  the  number  of  bilayers  N  increases  while  the  satiuration  field  resistivity 
Psat  remains  roughly  constant.  In  (b),  the  CPP-GMR  increases  from  about  25% 
to  40%  when  N  increases  from  about  17  to  40.  This  is  a  60%  change  in  the  GMR. 
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clearly  seen  from  the  difference  in  resistivities  which  shows  up  in  the  GMR  data  in 
4.2(b).  If  one  concludes  from  the  Figure  4.2(a)  that  Psat  is  roughly  constant,  then 
Po  must  be  an  increasing  fimction  of  A^.  By  detail  analysis  on  the  property  of  the 
samples,  the  authors  show  that  the  long  length  scale  (about  10  nm)  roughness  of 
the  interfaces  increases  as  A'^  increases,  and  suggest  that  the  increase  in  po  but  not 
Psat  is  due  to  the  increase  in  spin  dependent  scattering  at  the  interface.  Since  we 
only  include  atomic  scale  interface  roughness  in  our  calculation,  we  would  not  give 
the  same  trend  as  the  experiments. 

4.3    Comparing  Theory  with  Experiments 

In  Figure  4.4,  we  have  plotted  both  the  data  from  Gjis  et  al.  experiment 
and  those  from  the  Cyrille  et  al.  experiments  on  the  same  graph.  Cyrille's  data 
appears  to  be  on  the  y-axis  because  they  are  measured  at  low  temperatures,  which 
is  required  for  the  leads  to  be  superconducting.  As  we  can  see,  the  Gjis'  resistivity 
is  about  a  factor  of  two  larger  than  the  corresponding  data  of  Cyrille. 

Our  calculations  are  at  zero  temperature  so  we  do  not  explicitly  have  temper- 
ature effects;  however,  within  our  model,  we  have  several  disorder  parameters.  By 
increasing  the  bulk  disorder  parameter,  we  can  increase  scattering,  which  is  one  of 
the  effects  of  raising  the  temperatmre.  In  Figure  4.3,  we  have  plotted  the  the  max- 
imum and  saturation  resistivity  for  our  calculations  as  the  spin-independent  bulk 
disorder  varies.  A  moderate  residual  bulk  disorder  is  also  added  to  the  system. 
The  residual  bulk  disorder  is  chosen  so  that  resistivity  calculations  of  bulk  Fe  and 
Cr  films  agree  with  experiments.  The  results  in  Figiure  4.3  are  close  in  magnitude 
to  Cyrille's  data  and  about  twice  larger  than  Gijs'  data. 
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Figure  4.4:  Resistivity  in  Cyrille's  [Fe(30  A)/Cr(12  A)]n  samples  and  Gijs'  [Fe(30 
A)/Cr(10  A)]ioo  sample  plotted  on  the  same  graph.  The  saturation  field  (triangles) 
and  zero  field  (squares)  resistivity  in  Cyrille's  experiment  are  close  the  the  y-axis 
because  the  experiments  is  done  at  the  Uquid  helium  temperature.  The  saturation 
field  (pluses)  and  zero  field  (circles)  resistivity  in  Gijs'  sample  are  about  50% 
smaller  than  those  in  Cyrille's  experiments. 

In  order  to  fit  Gijs'  results,  we  have  to  calculate  the  resistivities  using  zero 
residual  bulk  disorder.  The  resistivities  calculated  this  way  are  close  to  Gijs'  data. 
However,  the  assmnption  that  the  bulk  of  the  sample  is  perfectly  grown  is  unrea- 
sonable. 

There  are  other  disorder  parameters  related  to  the  interface  that  we  can  change. 
In  this  calculation,  we  assumed  that  the  interface  thickness  to  be  4  monolayers, 
and  the  interface  mixing  percentage  to  be  20%.  The  resistivities  can  be  reduced 
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Figure  4.5:  Calculated  resistivity  in  [Fe(28.5  A)/Cr(11.5  A)]oo  for  moderate  resid- 
ual bulk  disorder.  Both  the  resistivity  for  parallel  (squares)  and  anti-parallel 
(crosses)  increase  with  the  bulk  disorder,  however,  the  parallel  increases  faster, 
therefore  the  GMR  will  decrease  with  the  bulk  disorder.  The  theory  data  is  close 
in  magnitude  to  Cyrille's  data  and  about  twice  larger  than  Gijs'  data. 


by  reducing  the  interface  disorder.  However,  we  doubt  that  the  interface  can  be 
much  sharper  than  4  monolayers  and  have  less  interface  mixing. 

A  better  way  to  look  at  the  differences  among  the  data  is  to  plot  the  change 
in  resistivity  against  the  saturation  resistivity,  as  shown  in  Figiure  4.3.  For  a 
fixed  temperature  and  sample,  there  are  two  independent  data,  the  resistivity  at 
zero  field  and  saturation  field.  By  plotting  functions  of  these  two  independent 
variables  against  each  other,  we  can  identify  trends  and  remove  the  uncertainty  in 
temperature  dependence  in  the  calculation.  As  seen  in  Figure  4.3,  Gijs'  data  and 
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the  theory  have  the  same  trend,  due  to  a  change  in  temperature  or  a  change  in 
bulk  disorder,  Cyrille's  data  represents  the  change  in  number  of  bilayers  and  has 
a  different  trend. 
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Figure  4.6:  Change  in  resistivity  against  saturation  resistivity.  The  data  of  Gijs 
(crosses)  and  the  theory  (squares)  have  the  same  trend;  while  the  data  of  Cyrille 
(circles)  has  a  different  trend  because  it  is  done  by  varying  the  number  of  bilayers. 


Another  advantage  of  plotting  the  data  this  way  is  when  we  compare  two  curves 
on  the  graph,  we  are  comparing  parameters  other  than  the  bulk  disorders,  which 
are  the  most  difficult  to  model  because  they  are  usually  not  directly  measured  in 
the  experiments.  With  the  most  uncertain  parameter  eliminated,  a  good  theory 
should  fit  the  experimental  ciurve  well. 
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As  seen  in  figures  4.4  and  4.3,  Gijs'  data  has  a  lower  resistivity  than  both  the 
theory  and  the  experiment.  The  trends  in  Figure  4.3  do  look  similar,  but  there 
is  a  large  offset  between  the  curves.  Fred  Sharifi  suggested  one  possibility  is  that 
there  may  be  other  conduction  paths  that  allow  electrons  to  by-pass  the  Fe/Cr 
superlattice.  In  Figure  4.3  we  can  obtain  Gijs'  curve  by  adding  a  shunt  resistor  to 
the  calculated  superlattice  model. 


Superlattice 


Figure  4.7:  Parallel  resistor  model.  The  effective  resistivity  is  found  by  1/peff  = 
1/Piatt  +  1/Pshunt5  where  Piatt  and  Pshunt  are  the  resistivity  of  the  superlattice  and 
the  shunt  resistor,  respectively. 

In  Figure  4.3,  we  plotted  the  the  change  in  resistivity,  po  —  Psat)  against  the 
saturation  resistivity  for  both  Gijs'  data  and  our  calculation  with  and  without  the 
shunt  resistor.  We  found  that  the  calculation  with  a  temperature-independent,  80 
ju!^cm  shunt  resistor  agrees  well  with  Gijs'  result.  Therefore,  Gijs  data  can  be  ex- 
plained if  in  the  experiment,  some  electrons  are  allowed  to  by-pass  the  superlattice 
and  thereby  reducing  the  measured  resistivity. 
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Figure  4.8:  Change  in  resistivity,  po  —  psat,  against  the  saturation  resistivity,  Psat- 
The  theory  data  without  the  shunt  resistor  (circles)  are  different  from  that  of  Gijs 
[Fe(30  A)/Cr(10  A)]  (crosses),  however,  the  theory  data  with  the  shunt  resistor 
(squares)  agrees  with  Gijs  results  very  well. 

In  order  to  fit  the  data,  the  resistivity  of  the  additional  path  has  to  be  almost 
independent  of  temperature  and  non-magnetic.  It  is  not  clear  what  caused  the 
extra  path.  One  possibility  is  that  pin-holes  are  formed  in  the  insulating  polyimide, 
which  allow  conduction  paths  of  Au  to  be  deposited  parallel  to  the  sample.  After 
the  polyimide  are  washed  away,  the  Au  conduction  paths  are  left  on  the  sample. 
Another  possibility  is  that  the  conduction  paths  may  come  from  impurities  at  the 
edge  of  the  superlattice. 

In  Figure  4.9,  we  compare  the  theory  with  Gijs'  [Fe(30A)/Cr(28A)]  data. 
Again,  with  the  shimt  resistor,  the  theory  agrees  with  the  experiment  very  well. 
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The  shunt  used  in  this  case  has  a  resistivity  of  75  /il2cm,  this  value  is  very  close  to 
the  one  obtained  for  the  Fe(30  A)/Cr(10  A)  sample.  This  suggests  that  the  two 
paths  are  made  of  similar  materials. 

By  using  our  superlattice  model  with  an  extra  conduction  path  way,  we  can 
explain  both  resistivity  and  GMR  data  for  two  different  samples  in  Gijs'  experi- 
ment. 
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Figure  4.9:  Change  in  resistivity,  po  -  Psat,  against  the  saturation  resistivity,  psat- 
The  theory  data  without  the  shunt  resistor  (circles)  are  different  from  that  of  Gijs 
[Fe(30  A)/Cr(28  A)]  data,  (crosses),  however,  the  theory  data  with  the  shunt 
resistor  (squares)  agrees  with  Gijs  results  very  well. 


81 


4.4    Other  Effects  on  the  GMR 

In  this  Chapter,  we  have  assumed  the  adjacent  magnetic  layers  to  be  completely 
anti-parallel  at  zero  magnetic  field.  In  experiments,  the  adjacent  magnetic  layers 
may  not  be  completely  anti-parallel  due  to  the  interlayer  coupling  [58-62].  The 
interlayer  coupling  controls  the  angle  between  the  orientation  of  magnetization  of 
adjacent  layers  at  zero  field.  In  Figure  4.10,  the  experiment  of  Parkin  shows  the 
saturation  field  oscillates  as  the  thickness  of  the  non-magnetic  spacer  layer  varies. 
This  effect  also  means  that  the  magnetic  coupling  oscillates  as  a  function  of  the 
layer  thickness.  The  origin  of  the  oscillation  is  the  RKKY  [58-59,  63-64]  interaction. 
The  GMR  is  also  shown  to  oscillate  with  the  thickness  in  the  experiment  of  the 
Michigan  State  University  group  [9-10].  The  effect  is  due  to  the  GMR  affected  by 
how  the  magnetic  layers  align. 


Spacer  Layer  Thickness  (A) 

Figure  4.10:  Dependence  of  saturation  field  on  spacer-layer  thickness  for  families 
of  Co/V,  Co/Mo,  and  Co/Rh  multilayers  [59]. 
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Figure  4.11:  The  magnetoresistance  of  Co/Cu  layers  as  a  function  of  the  thickness 
of  the  space  layer  (Cu)  for  current  perpendicular  to  the  plane  of  layers  (CPP) 
[9-10]. 


Another  complication  is  due  to  the  fact  that  the  magnetic  layers  are  adjacent 
to  non-magnetic  layers  or  the  vacuum.  In  such  case,  the  magnetization  in  both 
the  magnetic  layer  and  the  non-magnetic  layer  can  be  affected.  More  detailed 
information  can  be  reviewed  by  studying  the  spin-density,  such  as  those  by  Ohnishi 
et  al.  [42].  For  example,  Figm-e  4.12  shows  that  the  spin  density  on  the  Fe  layers 
extends  to  the  vacuum,  and  there  are  negative  spin  density  regions.  However, 
since  the  computational  cost  in  this  effect  is  high,  the  GMR  is  usually  computed 
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by  assuming  the  magnetic  layers  are  perfectly  aligned  in  the  saturation  field  and 
perfectly  anti-aligned  in  the  zero  field. 


SPIN  DENSITY  Fe  (OOl) 


Figure  4.12:  Self-consistent  spin-density  map  of  seven-layer  Fe(OOl)  in  unit  of 
0.0001  a.u.  on  the  (110)  plane.  Each  contour  line  differs  by  a  factor  of  2.  Dashed 
lines  indicate  negative  spin  density  [42]. 


CHAPTER  5 
OPTIMIZATION  OF  GMR  STRUCTURES 


In  this  chapter,  we  discuss  the  search  for  the  optimal  atomic  configuration  for 
the  giant  magnetoresistance.  In  previous  chapters,  we  have  discussed  how  to  model 
the  transport  properties  of  magnetic  nanostructures  and  compared  the  theory  to 
experiments.  Given  a  good  model,  we  should  be  able  to  make  predictions  as  well. 
In  particular,  what  atomic  configuration  gives  the  highest  GMR? 

This  question  is  difficult  because  the  configuration  space  is  large.  Take  mag- 
netic multi-layers  made  from  Sd  transition  metals  as  an  example,  suppose  we  were 
to  make  a  multilayer  that  consist  50  monolayers,  and  we  can  choose  to  use  among 
Co,  Ni,  or  Cu  atoms  in  each  monolayer,  there  are  roughly  7  x  10^^  possible  config- 
urations. It  is  hard  to  tell  which  configuration  gives  the  highest  GMR.  Although 
theories  are  capable  of  computing  the  GMR  for  any  of  the  configuration,  it  is  impos- 
sible to  exhaust  every  possibility.  The  conventional  way  of  solving  such  a  problem 
relies  on  physical  insights,  innovative  ideas,  and  often  trial-and-error  combined 
with  good  luck.  For  a  complicated  function  in  the  large  configm-ation  space,  such 
as  the  GMR,  explore  within  the  subspace  is  limited  to  one  that  is  simple  enough 
to  be  understood  with  physical  insights.  In  recent  years,  the  combinatorics  [65] 
method  has  been  a  useful  tool  in  searching  for  new  systems  for  interesting  physical 
properties.  For  example,  by  systematically  studying  samples  on  a  two-dimensional 
grid  of  gradual  gradient  of  compositions,  one  can  search  though  the  configuration 
space  at  a  much  faster  pace  than  was  previously  possible.  However,  trial-and-error 
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methods  are  unlikely  to  locate  solutions  in,  for  example,  the  50-dimensional  space 
mentioned. 

In  order  to  search  for  the  solution  in  a  high  dimensional  space,  the  com- 
puter program  must  not  search  the  configuration  space  randomly,  but  also  make 
use  of  the  local  information.  Optimization  methods  such  as  simulated  anneal- 
ing and  genetic  algorithms  can  be  used.  For  example,  using  simulated  annealing, 
Franceschetti  and  Zunger  [37]  are  able  to  find  the  largest  bandgap  configurations 
in  Alo.25Gao.75As  alloys  in  a  unit  cell  with  128  atoms  by  searching  for  about  10^ 
steps,  which  is  much  smaller  than  the  size  of  the  configuration  space. 

As  discussed  in  Chapters  2  and  3,  the  GMR  is  complicated  function  of  the 
atomic  configuration  and  the  disorder  parameters.  In  the  following  chapter,  we 
first  review  the  simulated  annealing  method.  Then  we  will  apply  the  method  to 
find  out  which  configuration  has  the  highest  GMR,  and  discuss  the  finding  and  its 
implication.  It  is  hoped  that  experimentally  one  can  grow  samples  according  to 
the  optimal  GMR  configuration  with  non-equilibrium  growth  techniques  such  as 
the  molecular  beam  epitaxy  (MBE). 

5.1    Inverse  Problems 

Let  /  be  a  function  of  the  atomic  configuration  C.  We  may  want  to  find  the 
minimum  (or  maximum)  of  the  function  in  the  configuration  space 


/rain  =  niin/(C), 


(5.1) 
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and  the  corresponding  configuration  Cmin-  In  other  situations,  we  may  want  to 
find  the  configuration  that  satisfies  the  target  condition  /(C)  =  /*^set^  which  can 
be  found  be  minimizing  of  the  absolute  value  of  difference.  More  complicated  re- 
quirements, such  as  prescribing  multiple  physical  property,  can  also  be  formulated 
as  the  minimization  of  the  function  [37] 

0{C)   =  E^.l/i(^)-/r''l'  (5-2) 
j 

where  Wj  is  the  weight  assigned  to  the  distance  of  property  fj  and  the  target  /j^"^^^*. 

This  problem  belongs  to  a  large  class  of  "inverse  problems" .  These  problems 
axe  difficult  because  the  functions  in  a  high  dimensional  configuration  space  are 
complicated,  and  intuitive  pictures  are  often  not  available.  Finding  the  highest 
GMR  configuration  is  an  example  of  an  inverse  problems.  In  order  to  solve  the 
inverse  problem  numerically,  one  must  evaluate  the  function,  or  solve  the  "direct 
problem,"  many  times.  In  practice,  the  direct  problem  itself  is  often  difficult  to 
solve.  For  example,  it  costs  about  5  x  10^^  floating  point  operations  to  find  the 
GMR  of  a  superlattice  with  a  unit  cell  of  50  atoms. 

To  solve  inverse  problems,  numerical  optimization  algorithms  such  as  simulated 
annealing  [32]  and  genetic  algorithms  [33-34]  are  used.  These  optimization  tech- 
niques have  been  used  to  solve  physics  problems  such  as  finding  the  equilibrium 
structure  by  minimizing  the  energy.  However,  the  optimization  describe  in  this 
chapter  is  rather  different  in  nature.  In  solving  the  inverse  problem,  the  solution 
obtained  satisfies  a  set  of  conditions  we  prescribe,  not  the  condition  set  by  natiire. 
This  allow  us  to  tailor-make  a  configuration  with  the  desired  properties  rather 


than  settle  for  what  is  naturally  available.  The  optimized  configuration  may  be 
realized  by  growth  techniques  that  can  build  structures  monolayer  by  monolayer, 
or  even  atom  by  atom.  In  addition,  the  configiiration  space  is  very  different  from 
those  used  in  the  finding  equilibrium  structmes.  Here,  we  can  change  the  chemical 
composition  in  the  configuration,  which  is  not  a  thermodynamical  process. 

5.2    Simulated  Annealing 

In  this  section,  we  explain  how  to  maximize  a  function  by  the  simulated  an- 
nealing algorithm.  The  flow  chart  of  the  algorithm  is  shown  in  Figiu-e  5.1.  To 
maximize  a  function  /  as  a  function  of  the  configuration  C,  an  initial  configura- 
tion C  is  randomly  generated  and  an  initial  "annealing  temperatiue"  T  is  set.  The 
value  of  f{C)  is  then  evaluated.  A  new  configuration  Ctriai  is  then  generated  by 
a  Monte  Carlo  move  from  the  previous  configmation,  and  the  value  of  /(Ctriai)  is 
also  evaluated.  If  the  functional  value  of  /  for  the  new  configuration  is  higher  than 
that  for  the  previous  configuration,  it  is  accepted.  On  the  other  hand,  if  the  new 
GMR  is  lower,  it  is  only  accepted  with  a  probabiUty  of  exp( (/(Ctriai)  -  /(C))/T). 
The  process  is  continued  as  the  annealing  temperature  is  gradually  lowered  until 
there  is  no  further  change  in  the  configuration  and  the  final  configuration  is  a 
near-global  maximum.  This  process  resembles  the  annealing  process  in  statistical 
mechanics,  where  the  total  energy  is  minimized  when  a  system  is  cooled  slowly. 

In  simulated  annealing,  the  value  of  /  is  used  to  decide  whether  to  keep  or 
discard  the  trial  configuration.  Compared  with  random  walks  routines,  simulated 
annealing  causes  the  system  to  spend  more  steps  near  the  configurations  possessing 
large  values  of  /.  Simulated  armeahng  is  easy  to  implement,  and  it  is  typical  that 
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the  number  of  configuration  tested  in  simulated  annealing  be  much  smaller  than 
the  size  of  the  configiuration  space. 
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Figure  5.1:  Flow  Chart  for  Simulated  Annealing.  The  flow  chart  for  using  simu- 
lated annealing  to  search  for  the  optimal  configm-ation  C,  that  maximize  a  function 
/,  is  shown. 


5.3    Accuracy  and  Speed  Requirements  on  Model 

Searching  for  the  optimal  structure  using  optimization  techniques  requires  a 
delicate  balance  between  acciuacy  and  speed.  On  one  hand,  the  model  has  to 
be  accurate  enough  so  that  the  optimal  solution  is  reliable.  The  model  has  to 
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capture  the  physics  and  to  accommodate  a  wide  range  of  configiirations.  The  level 
of  detail  in  the  model  varies  from  problem  to  problem.  For  optimization  of  the 
transport  properties  for  nanostructiires,  the  model  has  to  include  important  effects 
of  the  band  structure  of  the  materials,  the  geometry  of  the  nanostructures,  and  the 
scattering  effects  caused  by  impiirity  atoms  and  structural  disorders.  All  of  these 
factors  significantly  affect  electron  transport.  This  requirement  eliminates  methods 
containing  less  microscopic  detail  than  the  multi-band  tight-binding  model  that  we 
discussed  in  previous  chapters. 

On  the  other  hand,  calculations  of  the  transport  properties  in  this  model  have 
to  be  fast  enough  so  that  it  can  finish  in  reasonable  time.  The  number  of  anneal- 
ing steps  is  expected  to  be  of  the  order  of  10^  to  10^.  To  finish  the  optimization 
within  a  month,  one  has  to  be  able  to  compute  the  function  for  a  configuration 
in  about  4  to  40  minutes.  This  is  very  demanding  even  for  running  tight-binding 
codes  on  a  computer  cluster.  At  first  glance,  the  matrix  operations  of  finding 
the  Green's  function,  vertex  corrections,  and  conductivity  in  Chapter  2  are  are  all 
0{N^)  operations,  where  A'^  is  the  number  of  atoms  in  the  unit  cell.  In  addition, 
the  self-energy  and  the  vertex  corrections  are  solved  by  self-consistent  iterations. 
Performing  optimization  with  a  self-consistent  0{N^)  calculation  would  take  about 
50  CPU-years,  which  is  impractical.  This  is  why  we  have  developed  the  fast  algo- 
rithm in  Chapter  3.  For  A''  =  50,  this  algorithm  is  more  than  100  times  faster  than 
the  0{N^)  algorithm.  Running  on  a  Linux  cluster  with  ten  866MHz  Pentium  III 
CPUs,  it  takes  about  20  minutes  for  each  OMR  calculation  and  about  two  weeks 
for  an  optimization  with  1000  annealing  steps. 
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The  computational  cost  counting  in  Chapter  3  are  useful  in  estimating  the 
computer  resources  required.  A  fair  estimation  of  the  time  for  an  optimization  is 
required  before  coding  could  be  started. 

5.4    Constrained  Optimization  of  the  GMR 

In  this  section,  we  explain  how  to  apply  simulated  annealing  optimizing  the 
giant  magnetoresistance.  As  discussed  in  Chapter  3,  there  are  two  definition  of  the 
giant  magnetoresistance.  Here  we  adopt  the  "inflated"  definition  of  the  GMR: 

GMR  =   — -1,  (5.3) 

because  it  is  "optimization  friendly".  If  we  were  to  use  the  alternate  definition, 
GMR  =  1  —  aAp/crp,  the  GMR  would  be  bounded  by  100%  and  a  significant 
improvement  from,  say,  80%  and  85%  would  look  small  compared  to  the  aimeahng 
temperature,  and  could  be  washed  away  until  the  temperature  is  very  low,  hence 
extending  the  annealing  time  significantly.  When  the  physical  property  to  be 
optimized  is  "optimization  unfriendly,"  one  needs  to  find  a  monotonic  function  of 
the  property  that  magnifies  the  regime  of  interest. 

Not  every  configuration  qualifies  as  a  GMR  structure.  Unqualified  configura- 
tions, such  as  those  with  no  magnetic  layers,  will  have  zero  or  very  small  GMR. 
When  the  computational  cost  for  each  configuration  is  small,  one  can  rely  on  the 
exponentially  small  probability  to  make  sure  the  unquahfied  configurations  are  not 
chosen.  However,  when  the  computational  cost  for  each  configuration  is  large,  the 
time  wasted  in  computing  the  unqualified  structures  become  very  costly.  In  com- 
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puting  the  transport  properties,  self-consistencies  in  the  self-energy  and  the  vertex 
correction  are  required,  therefore  the  result  for  the  last  configuration  is  not  useful 
in  finding  the  result  for  the  present  configuration.  Therefore  each  GMR  calculation 
is  expensive.  One  can  also  use  a  test  algorithm  to  reject  unqualified  configurations, 
but  it  would  still  waste  significant  amount  of  time  when  the  qualifying  configura- 
tion subspace  is  much  smaller  than  the  full  configuration  space.  Therefore,  it  is 
better  to  perform  the  optimization  in  the  constrained  subspace  of  GMR  structures. 
The  algorithm  of  generating  the  configurations  in  the  constrained  subspace  will  be 
more  comphcated,  but  it  will  produce  results  in  a  much  shorter  time. 

To  optimize  the  GMR  using  simulated  annealing,  we  consider  the  constrained 
subspace  of  configurations  consisting  of  a  unit  cell  with  two  magnetic  layers  sep- 
arated by  two  non-magnetic  spacer  layers.  After  an  arbitrarily  chosen  initial  unit 
cell,  each  subsequent  configmation  is  generated  by  one  of  the  following  Monte 
Carlo  moves:  (i)  inserting  a  monolayer,  (ii)  removing  a  monolayer,  or  (iii)  chang- 
ing the  composition  of  a  monolayer.  To  make  sure  we  do  not  leave  the  constrained 
subspace  of  GMR,  a  move  is  rejected  unless  it  will  not  eliminate  a  layer,  or  turn  a 
magnetic  layer  into  a  non-magnetic  layer  by,  for  example,  changing  the  atoms  from 
Co  to  Cu.  To  prevent  complications  such  as  pin  holes  and  magnetic  dead  layers, 
we  further  limit  the  minimum  thickness  of  a  layer  to  two  mono-atomic  layers. 

If  the  GMR  of  the  new  configuration  is  higher  than  the  previous  configuration, 
then  it  is  accepted.  On  the  other  hand,  if  the  new  GMR  is  lower,  it  is  only 
accepted  with  a  probabihty  of  exp(AGMR/T),  where  AGMR  is  the  change  in  the 
GMR  and  T  is  the  simulated  annealing  temperature.  The  process  is  continued  as 
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the  annealing  temperature  is  gradually  lowered  until  there  is  no  fiirther  change  in 
the  configuration  and  a  final  near-global  maximum  in  the  GMR  is  attained. 

5.5    Optimal  GMR  Configuration  in  CoNiCu  superlattices 

In  the  following,  we  discuss  the  parameter  and  results  of  a  study  looking  for  the 
optimal  configuration  for  superlattices  made  with  Co,  Ni,  and  Cu,  in  conditions 
similar  to  GMR  experiments  at  the  room  temperature.  To  make  the  study  more 
manageable,  we  only  consider  the  fcc(lll)  lattices. 

As  mentioned  in  Chapter  2,  the  tight-binding  parameters  for  Equation  2.1  are 
obtained  from  fits  to  density-function  calculations  where  up  to  second-nearest- 
neighbor  hoping  energies  are  included.  In  order  to  simulate  experimental  condi- 
tions, we  fixed  the  bulk,  spin-independent  scattering  parameter  ^buik  at  0.2  eV^. 
Using  this  parameter,  the  resistivity  for  bulk  Cu  is  calculated  as  3.1  /iQcm,  about 
twice  that  of  the  room  temperature  resistivity  of  Cu.  Therefore,  we  are  mod- 
eUng  a  disordered  metal  system.  The  spin-dependent  impurity  potential  caused 
by  interface  diffiision  is  found  by  fitting  the  tight-binding  model  to  a  supercell 
density- functional  calculations  of  an  impurity  atom  in  the  host.  The  values  for  the 
impiurity  potentials  obtained  are  as  follows:  for  the  majority  spin,  Awco/Cu  =1-7 
eV,  AwNi/Cu  =0.56  eV,  Awco/Ni  =0.56  eV;  and  for  the  minority  spin,  Attco/Cu  =2.9 
eV,  AtiNi/cu  =1-05  eV,  Auco/m  =0.64  eV.  The  interface  mixing  percentage  is  fixed 
in  the  follow  way:  for  each  monolayer,  10%  of  the  atoms  diffuse  into  each  of  the 
two  monolayers  above  and  below.  Using  the  above  parameters,  the  GMR  for 
[C03/CU3/C03/CU3]  superlattice  is  36%.  This  value  is  similar  to  experiments.  The 
exact  value  for  the  mixing  percentage  in  the  experiment  is  imknown  because  it  is 
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very  difficult  to  measure.  However,  as  we  will  see,  the  optimal  atomic  configuration 
found  in  this  work  is  not  very  sensitive  to  the  errors  in  the  parameters  because  the 
effect  is  robust. 

Figure  5.2  shows  the  GMR  as  a  function  of  the  simulated  annealing  step.  The 
initial  configuration  (a)  is  [C03/CU3/C03/CU3],  a  symmetric  structure  made  of  four 
layers,  each  made  of  three  monolayers  of  either  Co  or  Cu.  Its  GMR  is  calculated  to 
be  36%.  However,  the  initial  configuration  does  not  affect  the  end  results.  Initially, 
the  annealing  temperatiue  T  is  set  to  be  60%.  As  the  optimization  progresses,  T 
decreases  slowly  and  reducing  the  probability  of  accepting  configiirations  which 
decrease  the  GMR.  As  we  can  see  from  Figure  5.2,  the  GMR  goes  up  to  about 
200%  and  then  goes  down  to  less  than  50%.  But  at  step  89,  the  GMR  jumps  to 
450%  and  no  other  configiirations  are  accepted  because  any  change  in  the  config- 
uration would  have  reduced  the  GMR  too  much  compared  with  T.  We  continued 
the  optimization  to  600  steps  and  no  fiurther  change  was  observed.  The  atomic 
configuration  for  the  highest  GMR  (d)  is  [Ni2/Cu2/Ni2/Cu2],  its  GMR  is  450%. 
This  is  a  very  surprising  result  because  most  room  temperature  studies  to  date 
have  GMR  less  than  100%.  The  GMR  is  sensitive  to  the  band  structure  and  hence 
the  atomic  configurations.  For  example,  configurations  (b)  is  [Ni2/Cu3/Ni2/Cu2], 
its  GMR  is  221%;  configuration  (c)  is  [Ni2/Cu2/Ni2/Coi/Cu2],  its  GMR  is  38%. 
Both  configurations  differ  from  the  optimal  configuration  (d)  by  only  one  mono- 
layer; however,  their  GMR  values  are  drastically  different. 
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Figure  5.2:  GMR  as  a  function  of  the  simulated  annealing  step.  An  initial 
configuration  (a),  which  does  not  affect  the  result,  is  first  chosen.  To  opti- 
mize the  GMR  with  simulated  annealing,  new  configurations  are  generated  by 
Monte  Carlo  moves  while  the  annealing  temperature  T  is  slowly  lowered.  When 
there  is  no  more  change  in  the  GMR,  the  process  stops  and  the  final  configura- 
tion (d)  is  the  optimal  configuration.  The  four  marked  configurations  and  their 
GMR  values  are:  (a)  [C03/CU3/C03/CU3]  36%,  (b)  [Nis/Cus/Nia/Cua]  221%,  (c) 
[Ni2/Cu2/Ni2/Coi/Cu2]  38%,  and  (d)  [Nis/Cus/Nia/Cua]  450%.  Notice  that  the 
GMR  is  very  sensitive  to  the  configuration.  Both  (b)  and  (c)  differ  from  (d)  by 
only  one  monolayer;  however,  their  GMR  values  are  very  diflFerent. 

5.6    Understanding  the  Large  GMR 


In  Chapter  4,  we  show  that  impurity  atoms  and  structural  disorder  often  reduce 
the  giant  magnetoresistance  significantly.  Therefore,  it  is  surprising  to  have  large 
GMR  for  similar  the  disorder  parameters.  We  note  that  the  GMR  as  a  function 
in  the  high  dimension  configuration  space  is  complicated.  Therefore,  it  is  possi- 
ble that  some  regime  of  the  configuration  space  has  been  overlooked.  To  get  a 
better  picture  of  the  GMR  as  a  function  in  the  configuration  space,  we  study  the 
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statistical  distribution  of  the  GMR.  To  do  this,  we  compute  the  GMR  of  a  large 
number  of  randomly  generated  configurations  that  satisfies  the  criteria  used  in  the 
last  section.  In  order  to  obtain  a  reasonably  accurate  distribution,  the  number 
of  configuration  in  a  statistical  study  is  a  lot  larger  than  the  number  of  steps  in 
the  simulated  annealing.  We  study  two  sets  of  configurations,  each  has  300  con- 
figiurations,  their  statistical  distributions  are  similax.  Therefore,  we  believe  the 
number  of  configurations  is  large  enough  to  give  an  reliable  pictmre  of  the  GMR 
distribution.  On  our  Linux  cluster,  it  takes  about  two  weeks  for  the  calculations. 

The  histogram  of  the  GMR  of  the  600  configurations  is  plotted  in  Figxure  5.3. 
The  average  GMR  is  22%,  and  the  standard  deviation  is  36%.  The  most  important 
feature  of  this  histogram  is  that  there  are  two  groups  of  configurations.  About  97% 
of  the  configmrations  are  located  near  the  main  peak,  they  have  GMR  less  than 
100%.  On  the  other  hand,  about  3%  of  the  configurations  have  GMR  many  stan- 
dard deviations  away  from  the  main  peak.  They  scatter  in  the  range  from  100% 
to  450%.^  Most  of  the  large  GMR  configurations  are  ultra-thin  Ni/Cu  superlat- 
tices,  only  two  of  the  large  GMR  configurations  contain  Co.  Since  the  large  GMR 
configurations  are  rare,  they  may  have  been  overlooked  in  previous  studies. 

The  very  different  statistical  distributions  in  the  two  group  suggests  there  is  a 
significant  physical  reason  behind.  In  order  to  understand  the  difference,  in  Figure 
5.4  we  study  two  typical  samples  in  each  group,  [Ni2/Cu2/Ni2/Cu2]  (circles)  for 
large  GMR,  [C02/CU2/C02/CU2]  (crosses)  for  regular  GMR.  We  plotted  the  Figure 
5.4(a)  £?-density  of  states  for  spin-up  channel  of  the  parallel  magnetic  configiuration, 

^The  optimal  configuration,  which  was  found  in  the  last  section  and  has  a  GMR  of  450%,  was 
not  among  the  600  randomly  generated  configurations  because  of  the  large  munber  of  possible 
configurations. 
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Figure  5.3:  Histogram  of  the  GMR  in  the  configuration  space.  In  order  to  get  a 
better  picture  of  the  GMR  in  the  configuration  space,  we  study  the  statistics  of  the 
GMR  for  600  randomly  generated  configurations.  The  average  GMR  is  22%,  and 
the  standard  deviation  is  36%.  As  we  can  seen,  about  97%  of  the  configurations 
are  located  near  the  main  peak.  They  have  GMR  less  than  100%.  On  the  other 
hand,  the  GMR  of  about  3%  of  the  configurations  are  many  standard  deviations 
away  from  the  main  peak,  they  scatter  in  the  range  from  100%  to  450%.  Most  of 
the  large  GMR  configurations  are  ultra-thin  Ni/Cu  superlattices.  Only  two  of  the 
large  GMR  configurations  contain  of  Co.  Since  the  large  GMR  configmrations  are 
rare,  they  may  have  been  overlooked  in  previous  studies. 

Figure  5.4(b)  the  conductivity  for  the  same  channel,  and  Figiu-e  5.4(c)  the  GMR, 
as  fmictions  of  the  Fermi  energy  in  the  rigid  band  picture.  The  physical  Fermi 
energy  corresponds  to  £'f  =  0  Ry.  As  seen  in  Figure  5.4(a),  the  two  configurations 
have  similar  band  structures.  The  main  difference  is  that  while  [Ni2/Cu2/Ni2/Cu2] 
has  a  sharp  d-DOS  despite  of  the  10%  interface  mixing  and  0.2  eV^  bulk  scattering, 
[C02/CU2/C02/CU2]  has  a  spread  d-DOS  even  at  5%  interface  mixing  and  the  0.2 
eV^  bulk  scattering.  This  is  because  the  scattering  between  spin-up  Co  and  Cu 
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is  about  a  three  times  the  scattering  between  spin-up  Ni  and  Cu.  Therefore,  the 
scattering  in  Co/Cu  interface  per  atom  is  about  9  times  as  effective  as  the  scattering 
in  the  Ni/Cu.  As  a  result,  [C02/CU2/C02/CU2]  has  a  spread  d-peak  and  its  (i-DOS 
at  the  Fermi  energy  is  large,  while  [Ni2/Cu2/Ni2/Cu2]  has  a  sharp  c?-peak  and  its 
d-DOS  at  the  Fermi  energy  is  small.  As  seen  from  figiures  5.4(a)  and  5.4(b),  the 
conductivity  increases  as  the  c?-DOS  decreases.  This  is  because  when  d-DOS  is 
large,  s-d  hybridization  reduces  the  average  Fermi  velocity  significantly,  and  hence 
reducing  the  conductivity.  Since  the  conductivities  in  the  spin-down  channel  and 
in  the  anti-parallel  magnetic  configiu-ation  do  not  vary  as  much  as  this  channel, 
the  GMR  curve  roughly  follows  the  conductivity  curve  for  the  spin-up  channel 
of  the  parallel  magnetic  configuration.  We  notice  that  both  examples  eventually 
have  large  GMR  at  high  enough  energy,  however,  only  the  GMR  at  the  Fermi 
energy,  u;  =  0,  is  measiured  in  many  experiments.  As  a  result,  there  is  a  qualitative 
difference  in  the  GMR  of  [Ni2/Cu2/Ni2/Cu2]  and  [C02/CU2/C02/CU2]  despite  of 
their  similar  band  structure. 

^  5.7  Conclusion 

Using  a  highly  optimized  GMR  code,  we  have  used  the  simulated  annealing 
algorithm  to  search  for  the  optimal  atomic  configuration  for  CPP-GMR.  For  condi- 
tions similar  to  experiments  at  the  room  temperatmre,  we  foimd  that  [Ni2/Cu2/Ni2/Cu2] 
has  the  highest  GMR,  which  is  450%.  This  is  smrprising  because  GMR  measure- 
ments for  most  configurations  are  usually  much  lower  when  scattering  is  large.  We 
have  also  studied  the  statistical  distribution  of  GMR,  which  has  given  us  insight 
about  GMR  as  a  fimction  in  the  configuration  space.  The  configurations  are  found 
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Figure  5.4:  Comparison  of  [Ni2/Cu2/Ni2/Cu2]  (circles)  and  [C02/CU2/C02/CU2] 
(crosses)  superlattices.  Although  their  band  structure  is  similar,  the  NiCu  super- 
lattice  has  a  large  GMR  (450%),  while  the  CoCu  superlattice  has  a  smaller  GMR 
(50%).  This  difference  is  caused  by  the  spreading  of  the  rf-peak  from  lower  ener- 
gies to  the  Fermi  level,  u;  =  0.  (a)  The  d-density  of  states  of  the  spin-up  channel 
for  the  parallel  magnetic  configuration,  c?-DOS||,  of  the  [Ni2/Cu2/Ni2/Cu2]  at  the 
Fermi  level  is  much  lower  then  that  of  the  [C02/CU2/C02/CU2].  (b)  This  leads  to 
a  higher  conductivity,  a||,  for  the  same  spin  channel  of  [Ni2/Cu2/Ni2/Cu2]  than 
in  [C02/CU2/C02/CU2].  Since  the  conductivities  a^,  a-^,  cr||  are  relatively  fiat 
compared  with  an  near  the  Fermi  energy,  the  GMR  (c)  resembles  cr-fi  (b). 


to  be  in  two  groups,  97%  of  the  configurations  have  regular  GMR,  and  3%  of  the 
configurations  have  large  GMR.  Since  the  large  GMR  configurations  are  rare,  they 
might  have  been  overlooked.  We  also  show  that  the  large  GMR  is  caused  by  the 
relatively  sharp  d-peak  and  hence  low  o?-DOS  at  the  Fermi  energy  for  the  spin-up 
channel  of  the  parallel  configuration,  which  in  turn  causes  high  conductivity  in  the 
same  channel  and  a  large  GMR. 
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We  have  studied  more  than  a  thousand  configurations,  however,  the  cost  of  hu- 
man resource  is  low.  Once  we  learned  that  large  GMR  exists  in  a  certain  regime, 
it  took  us  Uttle  time  to  study  this  regime  and  find  out  the  physics  behind.  In  a 
sense,  the  computer  decides  what  to  look  for,  in  a  tedious  but  slightly  'intelligent' 
way,  yet  it  is  effective  in  terms  of  actual  cost  for  the  research  project.  It  would  be 
a  waste,  if  at  all  possible,  to  have  people  performing  a  search  the  way  the  com- 
puter did.  To  conclude,  the  optimization  algorithms  has  made  our  research  more 
effective.  The  effect  found  here  is  a  robust  effect  with  simple  physical  explanation, 
yet  we  would  not  have  found  it  with  conventional  discovery  process. 


CHAPTER  6 
CONCLUSION 


In  this  final  chapter,  we  present  a  brief  review  of  the  dissertation  and  discuss 
its  implications  as  well  as  future  directions.  In  Chapter  2,  we  explained  how  to 
obtain  the  parameters  for  the  multi-band  tight-binding  Hamiltonian  and  for  the 
disorder.  We  also  discussed  how  to  calculate  the  conductivity  within  the  model, 
including  the  vertex  correction  necessary  for  current  conservation.  In  Chapter  3, 
we  developed  a  highly  optimized  parallel  code  capable  of  finding  the  self-energy 
in  0{N)  and  the  conductivity  in  0{N'^).  In  Chapter  4,  we  discussed  the  some  of 
the  experiments  in  the  current-perpendicular-to-plane  giant  magnetoresistance  and 
applied  oiu:  method  to  configurations  similar  to  the  experiments.  Our  calculations 
agrees  with  the  resistivity  and  giant  magnetoresistance  data  of  Cyrille  et  al,  and 
we  can  explain  quantitatively  data  of  Gijs  et  al.  by  a  shunt  resistor  model. 

Although  we  only  studied  the  magnetic  multilayer  superlattices,  our  method 
can  be  applied  in  a  wide  variety  of  systems  and  more  complicated  geometries.  With 
small  modification,  this  code  should  also  work  for  non-periodic  systems,  such  as 
systems  with  conducting  leads  or  systems  with  complicated  interface  structures. 
Most  importantly,  the  size  dependence  will  remain  the  same.  Also,  the  fast  con- 
ductivity code  may  also  be  used  to  study  large  systems. 

In  Chapter  5,  we  turned  to  a  new  direction  in  our  research.  Using  the  simulated 
annealing  method,  we  searched  for  the  optimal  giant  magnetoresistance  structiues 
made  with  Co,  Ni,  and  Cu  monolayers.  We  found  that  some  ultra-thin  Ni/Cu  su- 
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perlattices  have  a  high  GMR,  as  much  as  450%,  even  when  the  bulk  and  interface 
disorder  are  moderate.  This  is  a  surprising  discovery  because  it  is  expected  that 
disorder  at  this  level  should  reduce  the  GMR  significantly.  The  large  GMR  can  be 
understand  by  noting  that  in  a  small  number  of  ultra-thin  Ni/Cu  superlattice,  it 
is  possible  to  have  a  low  majority  d-density  of  states  at  the  Fermi  level  even  with 
moderate  disorder.  In  3d  transition  metals,  low  d-density  of  states  at  the  Fermi 
level  implies  low  resistivity  in  this  channel,  and  hence  a  large  giant  magnetore- 
sistance.  The  effect  found  in  this  study  is  not  very  sensitive  to  the  parameters 
because  it  is  robust. 

Although  this  effect  is  large  and  has  a  clear  physical  explanation  behind,  it  was 
not  discovered  before  because  the  regime  of  this  effect  is  very  small  compared  with 
the  size  of  the  configuration  space.  With  the  simulated  annealing,  the  computer 
is  "intelligent"  enough  to  avoid  trapping  by  local  optima  while  using  the  local 
information  to  find  the  global  optima.  This  method  makes  our  research  more 
efficient  and  eventually  leads  to  findings  that  are  not  likely  to  be  discovered  with 
conventional  methods.  As  far  as  we  know,  this  is  the  first  optimization  study  of 
transport  property  at  the  atomic  level. 

By  shifting  some  tedious  task  to  the  fast  performing  computers  of  today,  one 
can  make  a  qualitative  difference  in  the  process  of  scientific  discovery.  This  is  an 
important  trend,  although  what  we  used  is  very  far  away  from  the  "engines  of 
design"  predicted  in  Eric  Drexler's  Engines  of  Creation  [66],  which  "will  be  able 
to  create  bold  new  design  without  human  help".  We  still  need  to  ask  the  right 
question,  formulate  the  physical  model,  verify,  understand,  and  make  use  of  the 
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results.  Nevertheless,  a  significant  amount  of  hiunan  resources  which  would  have 
spent  in  trial-and-error  is  reduced. 

The  program  used  here  may  also  be  apphed  to  optimize  the  GMR  for  hot  elec- 
trons, where  the  measurement  is  made  at  energies  higher  than  the  Fermi  energy. 
This  approach  may  also  be  applied  to  the  area  of  spin  transport  across  semicon- 
ducting/metaUic  interfaces,  where  the  band  mismatch  is  very  large  and  optimal 
spin  transport  structmes  may  not  be  find  intuitively.  With  a  suitable  model  Hamil- 
tonian  for  the  semiconductors,  it  may  be  possible  to  optimize  spin  transports  across 
semiconducting/metaUic  interfaces. 

We  expect  more  types  of  problems  can  be  studied  with  this  approach  as  faster 
computers  are  available. 


APPENDIX  A 
FAST  INVERSE  ALGORITHMS 


To  compute  the  Green's  function  efficiently,  it  is  essential  to  have  fast  inverse 
algorithms  that  converge  well.  If  the  size  of  a  tight-binding  Hamiltonian  is  A'^, 
we  would  like  to  make  use  of  the  sparsity  of  the  obtain  the  full  Green's  function 
in  0{N'^).  In  addition,  since  only  the  onsite  Green's  function  are  needed  to  find 
the  self-energy,  it  is  desirable  to  have  an  algorithm  that  gives  the  diagonal  of  the 
inverse  in  0{N).  The  matrices  of  interest  are  semi-hermitian^  tridiagonal  block 
matrices,  with  or  without  periodic  boundary  conditions.  Although  the  inverse 
of  such  matrices  is  common  in  physics  applications,  there  are  not  many  inverse 
algorithms  or  codes  targeting  this  type  of  matrix,  especially  with  periodic  boundary 
conditions. 

In  Section  A.l,  we  review  the  algorithm  discribe  by  E.  M.  Godfrin  in  1991  to 
find  the  inverse  of  block  tridiagonal  matrices,  without  periodic  boundary  conditions. 
It  costs^  0{N'^)  to  compute  the  full  inverse  and  0{N)  to  compute  the  diagonal  of 
the  inverse. 

^Matrix  M  is  a  semi-hermitian  matrix  if  M{j,j')  =  M{j',j)*  for  all  j  ^  j'.  For  example, 
M  =  u;  —  H  +  iri  is  a,  semi-hermitian  matrix. 

^Since  the  elements  in  block  matrices  are  themselves  matrices,  element  additions  are  much 
cheaper  than  multiplications,  and  are  hence  neglected  in  cost  counting.  This  method  of  counting 
operations  is  different  from  that  of  general  matrices,  where  both  a  multiplication  and  an  addition 
are  count  as  floating  point  operations(flops)  and  cost  the  same  in  many  computers  with  math- 
coprocessors.  In  this  chapter,  the  elements  are  m  x  m  matrices,  and  the  operation  cost  counts 
are  therefore  in  unit  of  flops. 
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In  Section  A. 2,  we  develop  a  method  to  solve  the  inverse  of  tridiagonal  block 
matrices  with  periodic  boundary  conditions.  In  many  inverse  algorithms  such  as 
the  Godfrin  algorithm,  it  is  crucial  to  have  fixed  boundaries  from  which  iterations 
start.  Although  the  periodic  boundary  condition  only  adds  two  blocks  to  the  orig- 
inal matrix,  it  makes  finding  the  inverse  much  more  difficult.  To  overcome  this 
difficulty,  we  use  the  Woodbury  formula  to  compute  the  change  in  the  inverse  of  a 
matrix  after  a  rank  one  change.  Using  both  Godfrin  and  Woodbury  methods,  we 
are  able  to  compute  the  inverse  of  a  block  tridiagonal  matrix  with  periodic  bound- 
ary condition  and  the  diagonal  of  the  inverse  in  0(N'^)  and  0{N),  respectively. 

In  Section  A. 3,  we  improved  Godfrin's  method  for  inverse  for  both  types  of 
matrices.  Detailed  operation  counts  are  provided  to  show  the  importance  of  the 
results  presented  in  this  chapter. 

The  methods  in  this  appendix  are  tested  against  Matlab  for  large  matrix  in- 
verses to  double  precision  accuracy.  For  the  type  of  matrices  discussed  here,  these 
algorithms  should  converge  better  than  the  algorithms  for  general  matrices. 

A.l    Godfrin  Inverse  Algorithm  for  Tridiagonal  Block  Matrix 

Godfrin's  method  [67]  is  similar  the  one  derived  by  Economou-Cohen,  but  God- 
frin's result  is  exact  and  converges  very  well,  and  the  formulation  is  particularly 
suitable  for  our  calculations.  The  original  formulation  is  for  a  semi-hermitian  ma- 
trix, however,  it  can  be  formulated  for  a  general  tridiagonal  block  matrix. 

Let  M  be  a  block  tridiagonal  matrix,  and  G  =  be  the  inverse.  It  is  best 
to  demostrate  how  to  compute  the  inverse  with  a  4  x  4  matrix,  and  then  present 
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the  general  algorithm.  The  explicitly  definition  of  the  inverse  is  given  by 
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This  is  a  set  of  16  equations;  however,  the  equations  for  each  column  of  G  can  be 
solved  separately.  For  example,  the  equations  for  column  2  are: 


M11G12  +  M12G22  =  0  (A.l) 

M21G12  +  M22G22  +  M23G32  =  1  (A.2) 

M32G22  +  M33G32  +  M34G42  =  0  (A.3) 

M43G32  +  M44G42  =  0.  (A.4) 


This  set  of  equations  can  be  solved  independent  of  other  equations.  From 
Equation  A.4, 

G42     =  -M4^^M43G32 

M34G42   =  -M34M4-/M43G32 

=   -X3G32.  (A.5) 
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Substituting  the  result  into  Equation  A. 3,  we  have 
[Ms3  —  XzlGsi   =   —  M32G22 

^23^32     =     -M23[M33  -X3]-^M32G'22 

=   -X2G22.  (A.6) 

We  leave  behind  Equation  A. 2,  which  is  non-zero  on  the  right  hand  side,  until 
other  equations  are  used.  The  idea  is  that  we  first  represent  every  variables  above 
and  below  G22  in  terms  G22,  then  Equation  A. 2  can  be  solved  easily. 
Prom  Equation  A.l,  we  have 

G12   =   -Mfi^Mi2G22  .     ■  .  . 

M21G12   =  -M2iMf/Mi2Gi2 

=  -I2G22.     ^       :  /  , '  (A.7) 

Substitute  the  Equations  A.6  and  A. 7  into  Equation  A. 2,  we  obtain 

—  ^2^22  +  M22G22  —  X2G22     =     1  ' 

G22   =   [M22  -  X2  -  Fa]"'-  (A.8) 

The  rest  of  the  variables  are: 

G42  =  -M^M4iGz2 
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Gn   =   -Mf/Mi2Gi2.  (A.9) 

We  have  spent  0{N)  operations  in  order  to  solve  this  set  of  A'^  equations. 
However,  since  the  form  of  the  equations  on  the  left  side  of  Equations  A.1-A.4  are 
the  same  for  other  columns,  a  large  part  of  the  result  can  be  reused  when  solving 
equations  for  column  1,  3  and  4.  The  formal  procediu-e  is  as  follows  : 

Define  : 

Xn  =  0 

Xi  =  Mi,i+i[Mj+i,i+i  -  Xi+i]-^Mi+i^i  for  1  <  i  <     -  1 

n  =  0 

Yi  =  Mi,i_i[Mi_i,i_i  -  Fi_i]-iMi_i,i  for  2  <  i  < 

Ci  =  -[Ma  -  Xi]-'Mi^i_i  for  2  <  i  < 

A  =  -[Mii  -  Fi]-'Mi,i+i  for  1  <  i  <  A^  -  1,  (A.IO) 


then 


Gii  =  [Mu-Xi-Yi]-' 

Gij   =   CiGi-ij  for  i  >  j 

Gij   =   Did+ij.  for  i<j  (A.ll) 


It  only  takes  lOA^'  -  2  operations  to  find  X,  Y,  C,  D.  In  section  A.3,  we  will  see 
that  the  cost  of  finding  X,  F,  C,  D  can  be  reduced  to  6A^  -  2  by  rearranging  terms. 
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The  diagonal  and  any  column  j  of  Gij  can  be  found  in  0{N).^  The  cost  of  the 
full  inverse  is  A^^  +  0{N). 

A. 2    Inverse  of  Tridiagonal  Block  Matrix  with  Periodic  Boundary 

Condition 

Computing  the  inverse  of  tridiagonal  block  matrix  with  periodic  boundary  con- 
dition is  more  complicated.  We  are  still  able  to  find  the  inverse  in  0{N^)  and  the 
diagonal  of  the  inverse  in  0{N).  However,  the  coefficients  will  be  larger. 

In  the  following,  we  first  discuss  the  Woodbury  formula,  which  is  used  to  com- 
pute the  change  to  the  inverse  of  a  matrix  after  a  rank  one  change. 

A. 2.1    Rank-one  changes:  Woodbury  formula 

Let  M,  B  he  N  X  N  block  matrices,  u,  v  be  size  N  column  block  vector,  the 
Woodbury  formula  [68]  ^  states  that 

if       M  =  B-uv'^, 
then     M-^  =  B-^  +  B-\av'^B-\  (A.12) 

where  a  =  (1  —  v'^B~^u)~^  has  the  same  dimension  of  an  element  block.  To  prove 
this  formula,  one  considers  the  product 

{B-uv'^){B-^-^B-'uaii^B-^) 

^Note  that  in  order  to  find  a  row  in  0{N),  we  must  modify  the  algorithm,  as  discussed  in 
A.2.2. 

^It  is  called  the  Sherman- Morrison  formula  when  the  size  of  a  block  is  1  x  1. 
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uv 


''^B-^uav'^B- 


=   1  —  uv 


<'^B-^  +u{l-v 


<'^B-^u)av'^B 


-1 


% ... 


(A.13) 


because  {l-v'^B~^u)a  =  1.  Similarily,  the  product  {B~^+B~'^uav'^B''^){B-uv'^) 
is  also  equal  to  1.  Although  this  formula  looks  like  a  perturbation  method,  it  is 
indeed  exact. 

A.2.2    Applying  the  Woodbury  Formula 

The  Godfrin  algorithm  in  A.l  rehes  on  fixed  boundary  conditions  in  order  to 
start  the  iteration.  For  a.NxN  block  tridiagonal  matrix  M  with  periodic  boundary 
conditions,  there  are  two  more  block  elements:  Mijv  and  Mat i .  This  two  elements 
makes  the  inverse  much  harder  to  solve  than  it  looks  because  there  is  no  place  to 
starting  the  iteration. 

In  order  to  find  the  inverse,  we  first  write  M  into  a  tridiagonal  block  matrix 
without  the  boundary  condition  plus  the  periodic  conditions.  The  inverse  of  the 
first  term  is  computed  by  the  Godfrin  method,  and  the  correction  due  to  the 
periodic  condition  is  computed  by  the  Woodbm-y  formula.  In  other  words,  write 
M  in  the  following  form:  .  . 


M   =  B-uv'^ 
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Mu  -  MiN  Mi2   ...  0 
M21       M22  0 


0 
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1  0 


MiV-l,Ar-l  Mn~i,n 
Mn,n-i     Mnn  -  Mni  j 


0  1 


(A.14) 


where  B  is  the  tridiagonal  matrix  without  periodic  boundary  condition,  u  is  the 
column  vector,  and  is  the  row  vector.  B'^  can  be  found  by  Godfrin's  algorithm 
discussed  in  section  A.l.  can  then  be  computed  by  the  Woodbury  formula 

discussed  in  A.2.1. 

A. 2. 3    Cost  Considerations 

The  method  in  A. 2. 2  looks  very  simple,  however,  our  goal  of  computing  the 
diagonal  of  the  inverse  in  0{N)  can  be  acheived  only  by  examining  the  method 
carefully.  Since  u  and  v  in  Equation  A.14  are  sparse,  we  should  write  down  the 
Woodbury  formula  explicitly  for  our  case.  In  the  equation 


a  =  {l-vlBr^\)-\ 


(A.15) 


Ill 


where  repeated  indices  are  summed  automatically,  we  only  need  to  calculate 

-   {B^,^ +  B]:j\)ui  +  {B^j^UN  +  B]^l!)uN,  (A.16) 

It  costs  only  6  multiplications. 

The  correction  term  in  Equation  A.  10  is  more  complicated: 

{B-'uav''B-%  =  {B-'uaUv^B-%  (A.17) 

where 

{B~'^ua)i  =  Bii{uia)  +  BiNiuNOt), 
{v^B-%   =  vJ{B-%+vl{B-')^j 

=   {B-%  +  {B-')r,j.  (A.18) 

If  we  do  this  carefully,  i.e.  multiply  un  and  ui  by  q  first,  it  costs  AN  +  p  multi- 
plication, where  p  is  the  number  of  elements  we  want.  Therefore,  the  cost  for  the 
correction,  {B~^uav^B~^)ii,  to  the  diagonal  is  5N,  and  the  cost  for  the  correction 
to  the  whole  matrix  is  N'^  +  AN. 

In  Equation  A.16,  we  need  both  columns  {B''^)ii  and  (J5~^)ijv,  and  rows  {B~^)ij 
and  {B~^)f^ij.  However,  in  section  A.l,  we  only  know  how  to  obtain  {B~^)ij,  where 
j  varies,  in  a  row  form.  If  we  want  to  obtain  {B~^)ii^,  then  for  each  i,  it  takes 
0{N)  steps,  and  the  whole  column  costs  0{N'^).  In  order  to  obtain  the  columns 
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also  in  0{N),  we  need  to  find  the  inverse  of  the  tridiagnonal  matrix  in  the  cohimn 
form.  Similar  to  Equations  A.  10  and  A.ll,  we  have 


q  =  -Mi_,,i[Mu-X, 


1-1 


for  2<i<N 

for  l<i<N  -1,  (A.19) 


and, 


Gij  —  Gij-iCj 
Gij  =  Gij+iD'j 


for  i  <  j 
for  i  >  j. 


Note  that  Cj-  /  (Cj)^  because  Ma  is  usually  semi-hermitian  or  worse. 


(A.20) 


A. 3    Efficient  Coding  for  the  Godfrin  Algorithm 

We  can  speed  up  the  original  formula  listed  in  Godfrin's  paper  by  being  careful 
not  to  repeat  some  of  the  calculations.  Also,  both  the  column  and  the  row  Godfrin 
inverse  can  be  computed  at  the  same  time.  The  proper  procedure  is  as  follows: 


Xn  =  0, 

Wr,  =  [Xn-Mmn]-\ 

Cn  =  WnMn,n-u 

C'j^  =  Mn-i,nWn 
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and 

Y,   =  0, 

D,   =  Z,Mn, 

D[   =  M21Z1 

Y2   =  M2i[Mn-Yi]-'Mu 

=  -M21D, 

=  -D[Mn. 


(A.21) 


To  calculate  Xj,Yj,Wj,Zj,Cj,  Dj  for  all  j,  it  takes  6A''  —  2  multiplication  or 
inverse  operations.  It  takes  A'^  more  inverses  to  find  out  the  diagonals.  Therefore, 
the  cost  for  finding  the  inverse  of  a  tridiagonal  block  matrix  by  Godfrin  algorithm 
is  p  +  67V  +  0(1),  where  p  is  the  number  of  elements  we  need. 

It  takes  2N  more  calculations  to  find  all  the  C'^Dj,  and  4A'^  —  2  operations 
to  find  out  Gij,GNj,Gii,GiN-  Counting  the  cost  of  both  the  Godfrin  and  the 
Woodbury  methods,  the  total  cost  to  calculate  p  elements  of  the  inverse  of  a 
tridiagonal  block  matrix  with  periodic  boundary  condition  is  2p  +  16N  +  0(1). 
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This  method  is  very  efficient  in  finding  the  diagonals  when  TV  is  larger  than  about 
18. 

Since  these  methods  are  designed  specifically  to  solve  the  tridiagonal  block 
matrices,  it  should  be  more  efficient  and  accurate  than  Gaussian  elimination  with 
or  without  pivoting.  To  find  the  inverse  for  tridiagonal  block  matrix  with  periodic 
boimdary  condition  and  9x9  element  blocks,  it  takes  the  Matlab  sparse  matrix 
inverse  rountine  roughly  9.4A^^  x  9^  flops,  as  opposed  to  roughly  2N'^  x  9"^  flops  in 
this  method. 


APPENDIX  B 
COHERENT  POTENTIAL  APPROXIMATION 


The  coherent  potenetial  approximation  (CPA)  [26,  49-50]  is  an  average  method 
for  computing  the  density  of  states  in  substitution  alloys  within  the  single- impurity 
approximation.  Two  commonly  quoted  approaches  are  the  effective  medium  ap- 
proach and  the  locator  approach,  which  lead  to  the  same  set  of  equations. 

Consider  an  atom  of  type  j  in  an  effective  medium  with  potential  Ve-  In  real 
space  matrix  notation,  the  Green's  function  for  the  effective  medium  is 


where  we  have  written  the  Hamiltonian  as  a  sum  of  the  offsite  part  i/offsitej  and 
the  onsite  part  14  which  is  adjusted  later.  The  T-matrix  due  to  the  atom  is 


Go=[E  -  //offsite  -Ve  +  IT]] 


-1 


(B.l) 


(B.2) 


(B.3) 


where 


aj  ^  [7  -  {V,  -  K)Go] 


-1 


(B.4) 
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If  we  consider  only  a  single-site  impurity  ^  and  onsite  impurity  potentials,  then  Tj 
is  non-zero  only  at  a  single  site.  This  simplifies  the  problem  significantly.  Now, 
instead  of  solving  a  problem  for  the  entire  disordered  bulk  system,  we  only  have  to 
deal  with  a  single-site  T-matrix  (a  single  matrix  block/element  in  the  full  matrix). 
The  effect  on  the  whole  medium  will  be  carried  out  by  adjusting  the  effective 
potential  and  re-applying  Equation  B.l.  w''-^ 

The  CPA  requires  the  ensemble  average  scattering  to  be  zero  at  each  energy 
E,  i.e. 


This  condition  can  be  satisfied  by  modifying  Vg.  The  new  Ve  can  be  found  by 
rewriting  the  above  equation  as 


For  a  multiband  Hamiltonian,  the  order  of  operations  can  not  be  changed  because 
these  are  matrix  equations. 

Equations  B.l  and  B.6  are  then  applied  repeatedly  to  obtain  a  self- consistent 
solution  for  Ve  and  Gq.  The  above  algorithm  is  tested  to  work  for  a  9-band  tight- 
binding  model  of  fee  NiFe  and  bcc  FeNi.  The  effective  potential  K  looks  like  a 
weighted  average  of  Vi  and  V2  where  the  energy  dependent  weight  is  adjusted  by 
Go  and  K  itself.  Obviously,  the  above  equation  works  well  in  the  pure  host  limits, 
Ci  =  0  or  C2  =  0. 

^Thus  CPA  neglects  cluster  effects. 


(T)  =  ciTi  +  C2T2  =  c,ai{V,  -  K)  +  C2«2(^2  -  K)  =  0. 


(B.5) 


Ve  =  [Cifti  -I-  020:2]    ^{ciCliVi  +  0202^2). 


(B.6) 
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After  obtaining  the  effective  potential  and  Green's  function  of  the  effective 
medium,  we  can  compute  the  Green's  function  for  a  specific  sites.  For  an  atom  of 
type  j,  the  Green's  fimction  is 

Gj  =  Go  +  GoTjGo.  (B.7) 

This  equation  is  important  for  finding  the  site-specific  quantities  such  as  the  local 
density. 
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Magnetic  nanostructures  are  magnetic  devices  which  have  at  least  one  dimen- 
sion on  the  order  of  one  nanometer,  10"^  m,  and  they  have  become  the  focus  of 
a  new  trend  in  nanoelectronics  because  of  their  applications  in  magnetic  sensors, 
magnetic  memories,  and  magnetic  reading  heads  in  hard  disks  drives.  For  exam- 
ple, a  reading  head  in  a  magnetic  hard-disk  is  a  nanostructure  made  by  growing  a 
magnetic  layer  on  a  substrate,  followed  by  a  non-magnetic  layer,  and  then  another 
magnetic  layer.  The  thickness  of  these  layers  are  of  the  order  of  one  nanometer. 
The  resistance  of  this  structure  changes  as  a  magnetic  field  is  applied,  therefore  it 
is  used  as  in  a  reading  head. 

In  this  dissertation,  we  first  develop  a  method  to  model  the  electrical  response 
of  the  magnetic  nanostructure  in  the  atomic  scale.  This  method  is  then  improved 
such  that  it  can  be  applied  in  very  complicated  structures.  Then  we  apply  this 
method  to  study  the  giant  magnetoresistance  effect  and  successfully  explain  some  of 
the  existing  experiments.  Finally,  we  turn  to  a  new  direction  of  research  in  nanos- 
tructures in  which  we  use  an  optimization  method  on  high  performance  computers 


to  search  for  new  nanostructures.  We  found  a  magnetic  nanostructure  that  has  a 
surprisingly  large  giant  magnetoresistance  ratio.  It  is  hoped  that  this  nanostruc- 
ture will  be  realized  in  experiments.  Using  computers  to  search  for  thousands  of 
nanostructures  automatically,  the  efficiency  of  scientific  discovery  is  increased  and 
a  significant  amoimt  of  human  resources  is  saved. 
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